Análisis sistemático de variantes genéticas y epigenéticas funcionales en cáncer colorrectal.

El cáncer colorrectal (CCR) es una neoplasia genéticamente compleja impulsada por alteraciones genómicas y epigenómicas acumulativas. Si bien las mutaciones codificantes (por ejemplo, APC y KRAS) son factores desencadenantes bien caracterizados (^[1]), más del 90% de los variantes asociados a la enfermedad identificados en estudios de asociación del genoma completo (GWAS) se encuentran en regiones no codificantes de proteínas (^[2]), lo que sugiere funciones críticas para la desregulación de los elementos reguladores. Los potenciadores, en particular, coordinan la expresión génica espacio-temporal a través del bucle de la cromatina (^[3]), y su alteración por polimorfismos de un solo nucleótido (SNP) o metilación del ADN contribuye a la oncogénesis (^[4]).

Los potenciadores exhiben una profunda especificidad de tejido y tipo celular en sus patrones de actividad, una propiedad fundamental para la regulación génica espacio-temporal precisa durante el desarrollo y la enfermedad (^[3]). En el cáncer, esta especificidad impulsa programas oncogénicos dependientes del contexto tumoral, donde variantes genéticas idénticas pueden ejercer efectos divergentes en diferentes neoplasias (^[5]). Más del 70% de los variantes de riesgo de cáncer muestran efectos potenciadores específicos del tejido (^[6]). Por ejemplo, uno de los SNP de riesgo de CCR más estudiados, rs6983267, se mapeó a un superpotenciador en 8q24, un locus de "desierto génico" ubicado a unos 335 kb corriente arriba del promotor del protooncogén MYC (^[7]). Su alelo "G" reclutó la unión del factor de transcripción (TF) TCF7L2 para amplificar la oncogénesis impulsada por Wnt. Este estudio fundamental estableció el paradigma de los SNP no codificantes que regulan los protooncogenes a través de interacciones de largo alcance de la cromatina. Cabe destacar que, en el cáncer de páncreas y las células de leucemia mieloide aguda (LMA), los potenciadores funcionales se encuentran corriente abajo del gen MYC. Este hallazgo subraya la naturaleza específica del tejido de la función del potenciador y el impacto regulador de los variantes dentro de estos elementos (^[4]). Estudios recientes revelan que la progresión del CCR también implica una reprogramación epigenética dinámica de los potenciadores. La transición metastásica se caracteriza por una amplia usurpación de potenciadores (^[8]), mientras que la metilación aberrante en las islas CpG asociadas a los potenciadores silencia los supresores de tumores (^[9]). Los enfoques de genómica funcional, como el ensayo de reportero masivamente paralelo (MPRA) en células primarias (^[10]) y el cribado de variantes relacionadas con la actividad silenciadora (^[11]), han proporcionado información valiosa sobre la función molecular de los variantes genéticos no codificantes. Sin embargo, la anotación sistemática de la función reguladora de los variantes no codificantes y las regiones metiladas diferencialmente (DMR) en los potenciadores del CCR, especialmente en contextos metastásicos, sigue siendo limitada, ya que los métodos tradicionales, como los ensayos de reportero, generalmente carecen de escalabilidad.

En este estudio, primero realizamos un análisis pan-cáncer de las regiones de potenciadores activos en 20 neoplasias utilizando el ensayo de accesibilidad de la cromatina mediante transposición y secuenciación (ATAC-seq) y los datos de secuenciación de inmunoprecipitación de la cromatina con H3K27ac (ChIP-seq). Este marco comparativo permitió una discriminación sistemática de los variantes reguladores comunes con un impacto pan-cáncer de las alteraciones específicas del tumor únicas de la patogénesis del CCR, priorizando así los factores desencadenantes dependientes del contexto para el cribado funcional. Para los potenciadores del CCR, capturamos los SNP reguladores y los sitios de metilación CpG correspondientes. Posteriormente, diseñamos y sintetizamos bibliotecas de oligonucleótidos de alta complejidad que encerraban fragmentos genómicos de 120 pares de bases (pb) que rodeaban los SNP o los sitios CpG de interés. Aprovechando dos estrategias de alto rendimiento, la secuenciación de regiones reguladoras activas autorreplicantes de SNP (SNP-STARR-seq) (^[12], ^[13]) y la secuenciación de regiones reguladoras activas autorreplicantes de metilación (Methyl-STARR-seq) (^[14]), en células del CCR, incluidas dos líneas celulares epiteliales derivadas del CCR primario, HCT116 y SW480, así como una línea celular derivada del CCR metastásico, SW620, del mismo origen que SW480, cuantificamos sistemáticamente el impacto de los variantes genéticos y epigenéticos en la actividad del potenciador e identificamos 922 SNP comunes con alteraciones significativas en la actividad del potenciador compartidas entre las líneas celulares primarias (HCT116 y SW480), junto con 3136 SNP reguladores adquiridos en la metástasis únicos de SW620. De manera similar, Methyl-STARR-seq reveló 1763 y 3588 elementos potenciadores sensibles a la metilación en HCT116 y SW480, respectivamente, y 4112 en SW620, de los cuales 3008 se desregularon específicamente en el contexto metastásico. A través de la integración multiómica [Hi-C, locus cuantitativo de expresión (eQTL) y secuenciación de ARN (RNA-seq)], vinculamos estos variantes a sus genes diana reguladores.

Al dilucidar la función molecular de estos variantes, demostramos que el variante rs67941642 [un SNP en alto desequilibrio de enlace (LD) con el SNP líder de GWAS rs6061231 en el cromosoma 20q] podría mejorar de manera dominante la unión del represor transcripcional del factor de crecimiento independiente 1 (GFI1) y promover la expresión del gen supresor de tumores TAF4, cuya regulación a la baja impulsó la proliferación y la migración de las células del CCR. Para los variantes asociados a la metástasis, también revelamos que el alelo "A" del SNP rs1962004, que se adquirió somáticamente en la línea metastásica SW620, creó un sitio de unión para los componentes del activador de proteínas 1 (AP-1), los factores de transcripción Jun protooncogénico, AP-1 (JUN) y el protooncogénico JunD, el factor de transcripción AP-1 (JUND), activando el oncogén LRRC61 para promover los fenotipos metastásicos y correlacionándose con un mal pronóstico del paciente. En términos de los sitios de metilación relacionados con el riesgo de CCR, encontramos que la hipermetilación en cg08640619 interrumpió la unión del TF RUNX2, lo que condujo a la regulación a la baja de los genes cercanos KIRREL1 y ETV3. En esencia, tanto cg08640619 demostró un rendimiento excepcional como biomarcadores de detección temprana [área bajo la curva de características operativas del receptor (AUROC) > 0,96]. Nuestro estudio introduce una potente línea experimental para el cribado funcional de los variantes genéticos y epigenéticos reguladores.

RESULTADOS

Diseño del proyecto para identificar variantes funcionales en el CCR primario y metastásico

Cómo los variantes reguladores dependientes del contexto impulsan la patogénesis del CCR, particularmente durante la metástasis, no se había capturado adecuadamente en estudios anteriores. Aquí, desarrollamos una estrategia que pudo evaluar sistemáticamente tanto los variantes genéticos como los epigenéticos para su impacto en la actividad reguladora transcripcional de los elementos potenciadores del CCR (Fig. 1).

Primero, definimos un conjunto de regiones potenciadoras del CCR putativas identificando 18.880 intervalos genómicos que se superponían entre los picos de ChIP-seq de H3K27ac y ATAC-seq en la línea celular del CCR HCT116, lo que indica tanto la actividad reguladora activa como la accesibilidad de la cromatina. Para distinguir el paisaje regulador del CCR de los elementos pan-cáncer comunes, realizamos un análisis epigenómico comparativo en 20 neoplasias (tablas S1 y S2), lo que nos permitió clasificar estos potenciadores como comunes en muchos cánceres o más restringidos en su perfil de actividad.

Posteriormente, utilizamos dos sistemas distintos, SNP-STARR-seq y Methyl-STARR-seq, para evaluar la actividad del potenciador que contenía las variantes genéticas y epigenéticas, respectivamente. Para la biblioteca SNP-STARR-seq, sintetizamos una biblioteca genómica de 120 nucleótidos (nt) que abarcaba todos los sitios de SNP dentro de las 18.880 regiones potenciadoras putativas y los loci de SNP de GWAS y los SNP enlazados a GWAS. Para cada SNP, se incluyó una secuencia genómica de 120 pb centrada en el variante. Todos los posibles alelos variantes (A, T, C y G) en el sitio polimórfico se sintetizaron como oligonucleótidos individuales, lo que permitió una evaluación exhaustiva de los efectos específicos del alelo para los variantes genéticos humanos conocidos. Para la biblioteca de metilación, nos centramos en dos categorías de secuencias genómicas: (i) ventanas de 120 nt centradas en los picos de H3K27ac, ya que estas regiones de potenciador y promotor activas son puntos calientes establecidos para la reprogramación epigenética en el cáncer, donde los cambios en la metilación del ADN ejercen su impacto funcional más profundo. (ii) regiones genómicas correspondientes a los DMR del CCR. Tras la construcción de la biblioteca de plásmidos, realizamos un tratamiento enzimático in vitro utilizando la metiltransferasa M.SssI para generar dos bibliotecas distintas: biblioteca de plásmidos totalmente metilada (tratada con M.SssI) y biblioteca de plásmidos de control no metilada (tratada con un simulacro) (ver Materiales y Métodos). Las puntuaciones de actividad del potenciador (EAS) se derivaron de los recuentos de lectura de ARN/ADN normalizados por DESeq2 como EAS = (ARN + 1)/ADN. Los cambios de alelo o metilación se calcularon comparando los valores de EAS entre el alelo alternativo/de referencia o la metilación/no metilación, con significación estadística. A través del cribado sistemático en tres líneas celulares del CCR, modelos de adenocarcinoma primario (HCT116 y SW480) y una línea celular derivada de la metástasis (SW620), identificamos variantes reguladoras funcionales moduladas significativamente tanto por las mutaciones genéticas como por las alteraciones epigenéticas. Los SNP de actividad preferente (paSNP; definidos como SNP que confieren una puntuación de potenciador diferencial significativa y un valor de P ajustado < 0,05; ver Materiales y Métodos) y los sitios de metilación de actividad preferente (paMetilsitos; definidos como metilación CpG que confiere una puntuación de potenciador diferencial significativa y un valor de P ajustado < 0,05; ver Materiales y Métodos) identificados en las tres líneas celulares se presentan en las tablas S3 y S4.

Al integrar aún más los conjuntos de datos clínicos relevantes de RNA-seq, eQTL y epigenómicos, asignamos los variantes reguladores a sus genes diana. Centrándonos en los variantes implicados en el desarrollo y la metástasis del CCR, validamos la relevancia funcional de los variantes identificados obtenidos a través de esta línea de cribado utilizando técnicas que incluyen ensayos de reportero de genes, inmunoprecipitación de cromatina (ChIP) y edición del genoma basada en CRISPR. En conjunto, este marco de genómica funcional integrado no solo descifra los vínculos mecanicistas entre los variantes no codificantes y la patogénesis/metástasis del CCR, sino que también identifica objetivos terapéuticos previamente no identificados y biomarcadores de diagnóstico altamente sensibles con un potencial de traslación significativo.

SNP-STARR-seq identificó SNP que confieren actividades de potenciador diferenciales en las células del CCR

Dentro de los 18.880 potenciadores candidatos, seleccionamos los SNP con una frecuencia de alelo menor (MAF) ≥ 0,05 tanto en las poblaciones europeas como en las asiáticas orientales del Proyecto de los 1000 Genomas (^[15]), lo que resultó en 24.526 SNP comunes. Para investigar el potencial regulador de los SNP dentro de las regiones del potenciador, incorporamos datos del Catálogo de GWAS y recuperamos 238 loci asociados al CCR significativos a nivel del genoma, junto con 6125 SNP adicionales en LD [coeficiente de determinación (R2) ≥ 0,5] con estos loci. Al combinar los SNP identificados a partir de las regiones del potenciador del CCR con los SNP relacionados con GWAS, obtuvimos 30.790 SNP únicos como elementos reguladores candidatos para el cribado funcional. Para cada SNP, se extrajo una secuencia genómica de 120 pb centrada en el variante (± 60 pb). Todos los alelos variantes se sintetizaron como oligonucleótidos individuales. En total, diseñamos y sintetizamos 91.566 oligonucleótidos únicos, formando una biblioteca de SNP-STARR-seq de CCR para un ensayo de reportero masivamente paralelo (Fig. 2A). La agrupación jerárquica demostró una alta reproducibilidad en ambos, los replicados biológicos y los replicados técnicos [reacción en cadena de la polimerasa (PCR)] del mismo tipo de célula (Fig. 2B). Para cada SNP, el efecto específico del alelo sobre la actividad del potenciador se cuantificó calculando el cambio transformado logarítmicamente de EAS entre el alelo alternativo y el alelo de referencia. Todos los SNP calificados se clasificaron luego en función de este valor, y los SNP con actividades de potenciador diferencialmente significativas se seleccionaron para la caracterización funcional posterior (HCT116, n = 1969 SNP; y SW480, n = 1846 SNP; Fig. 2C y tabla S3). La validación de los resultados de SNP-STARR-seq utilizando un ensayo de reportero de luciferasa de un solo locus demostró una fuerte concordancia con los resultados de secuenciación de alto rendimiento en ambas líneas celulares (Fig. 2D y fig. S1A).

Tanto las líneas celulares HCT116 como SW480 representan adenocarcinomas colorrectales en etapa temprana, que exhibieron una superposición significativa de variantes funcionales (P = 0.048, prueba exacta de Fisher; Fig. 2E). A continuación, nos centramos en 922 SNPs comunes que mostraron efectos significativos en ambas líneas celulares (Tabla S5). La función de estas variantes no codificantes también está respaldada por datos de eQTL (82.65%), lo que demuestra su impacto regulador en la expresión génica proximal o distal (Fig. 2F). Además, se predijo computacionalmente que más del 75% de estas variantes albergan sitios de unión a factores de transcripción (Fig. 2G). Se obtuvieron resultados similares a partir de los análisis de cada línea celular individualmente (Fig. S1, B y C). El análisis de motivos mediante la optimización hipergeométrica de enriquecimiento de motivos (HOMER) reveló además un enriquecimiento significativo de factores de transcripción (FT) asociados con el CRC dentro de los segmentos genómicos que muestran diferencias significativas en la actividad del potenciador alélico, incluidos el factor de transcripción tipo ATF con cremallera leucina básica (BATF), FOS like 2, la subunidad del factor de transcripción AP-1 (FOSL2) y el antígeno relacionado con Fos-2 (FRA2), que se han relacionado directamente con la progresión y metástasis del CRC (^[16], ^[17]), junto con JUNB, un componente del complejo AP-1 que se sabe que está desregulado en el CRC (Fig. 2H y Fig. S1D) (^[18]). Basándonos en nuestra anotación de potenciadores pancancerosos mencionada anteriormente, el 81.6% de los 922 SNPs comunes residían dentro de las regiones de potenciadores putativos del CRC, aunque algunas variantes relacionadas con GWAS se encuentran fuera de estas regiones. Este enriquecimiento es significativamente mayor que el observado en otros 19 tipos de cáncer, lo que confirma que nuestro sistema de alto rendimiento identificó eficazmente variantes regulatorias relevantes para el contexto del CRC (Fig. 2I). En conjunto, estos hallazgos demostraron que el SNP-STARR-seq podría capturar con éxito la gran mayoría de las variantes funcionales.

El papel epistático del SNP rs67941624 en la modulación de la actividad del potenciador en las células de CRC

Entre los 922 SNPs putativamente funcionales, rs6061231 se encuentra en un locus de riesgo de CRC y se identificó como un SNP principal por GWAS (^[19]), mientras que el SNP rs67941642 está en completo desequilibrio de enlace (R2 = 1, https://ldlink.nih.gov/) con rs6061231 (Fig. 3A). Ambos SNPs residían dentro de las regiones genómicas que presentaban fuertes marcas epigenéticas asociadas con potenciadores (H3K4me1 y H3K27ac) en las células de CRC humanas (Fig. 3B). Dada la estrecha proximidad y el vínculo genético, utilizamos un ensayo de reportero de luciferasa dual para investigar el papel sinérgico de las dos variantes. Primero, confirmamos la actividad del potenciador de los fragmentos genómicos que contienen cualquiera de los alelos mayor o menor de los SNPs rs6061231 y rs67941642. En particular, el alelo rs67941642-T exhibió una actividad del potenciador significativamente mayor que el alelo C, y el alelo rs6061231-C mostró una actividad significativamente mayor que el alelo A (Fig. 3B), lo que es consistente con el resultado de SNP-STARR-seq (Tabla S5). En la cohorte combinada de Asia Oriental y Europa (EAS + EUR), el análisis de haplotipos del locus rs6061231/rs67941642 reveló frecuencias del 77.0% para el haplotipo C-C (rs6061231-C/rs67941642-C) y del 23.0% para el haplotipo A-T (rs6061231-A/rs67941642-T), mientras que los otros haplotipos (es decir, C-T y A-C) prácticamente no ocurren en la población. Dado que los dos SNPs en desequilibrio de enlace exhibieron actividades de potenciador opuestas, preguntamos qué SNP fue dominante en la determinación de la actividad del potenciador alélico. Por lo tanto, construimos los plásmidos de reportero que contienen los haplotipos C-C y A-T, respectivamente, y los transfectamos a las células de CRC para una comparación cuantitativa. El resultado del ensayo mostró que el haplotipo A-T exhibió una actividad del potenciador significativamente mayor en comparación con el haplotipo C-C (P < 0.001, Fig. 3C), lo que indicó que el SNP rs67941642 confirió un papel predominante en la modulación de la actividad del potenciador dentro de su contexto genómico. El SNP rs6061231 fue un SNP principal en GWAS, lo que sugiere que el SNP principal resultó ser un proxy de la variante funcional del SNP rs67941642. Si bien ambos SNPs exhibieron actividad del potenciador específica del alelo, el haplotipo A-T (rs6061231-A/rs67941642-T) confirió una actividad significativamente mayor que el haplotipo C-C (Fig. 3C). Esto indica que la producción regulatoria neta de este locus está dominada por el efecto combinado de ambos SNPs, con rs67941642-T desempeñando un papel predominante en la determinación de la fuerza general del potenciador de los haplotipos humanos comunes. Se ha establecido que algunos genes exhiben un papel epistático sobre otros genes en la expresión del fenotipo. Aquí, revelamos elementos reguladores epistáticos que mostraron un papel dominante sobre otros potenciadores en la regulación de la transcripción del gen diana. El resultado subrayó aún más la importancia de analizar exhaustivamente la función molecular de las variantes en desequilibrio de enlace y adyacentes a los resultados de GWAS.

Estudios previos han demostrado que numerosas variantes genéticas no codificantes funcionales pueden alterar la unión de los factores de transcripción. El análisis de motivos sugirió la unión específica de los factores de transcripción GFI1 a rs67941642-T, p53 a rs6061231-C y SIX homeobox 2 (SIX2) a rs6061231-A (Fig. 3D y Fig. S2A). La genotipificación confirmó la heterocigosidad de ambos SNPs rs67941642 y rs6061231 en las células HCT116 y SW480 (Fig. S2B). Posteriormente, la inmunoprecipitación de cromatina (ChIP) seguida de la reacción en cadena de la polimerasa cuantitativa (qPCR) y la secuenciación del amplicón ChIP reveló el reclutamiento preferencial de GFI1 a rs67941642-T (P < 0.001) y de p53 a rs6061231-C (P < 0.001), mientras que la unión de SIX2 no mostró un sesgo alélico significativo en rs6061231 (Fig. 3E y Fig. S2C). Dada los efectos contrastantes de ambos SNPs en la actividad del potenciador, investigamos si los dos factores de transcripción podrían tener un papel interferente en este locus. Por lo tanto, desactivamos individualmente p53 y GFI1 en las células HCT116 y encontramos que la unión de GFI1 en rs67941642 o de p53 en rs6061231 no se vio significativamente afectada, respectivamente, lo que sugiere que no existe una influencia mutua entre estos dos factores de transcripción (Fig. S2D).

El locus rs6061231/rs67941642 regula cis la expresión de TAF4

Un estudio previo demostró que la amplificación del número de copias del cromosoma 20q sirve como biomarcador pronóstico en el subtipo de CRC con microsatélites estables (^[20]). Dado que los SNPs rs6061231 y rs67941642 se encuentran dentro del locus Chr20q, nos propusimos comprender su relevancia para la significación clínica. Primero, pretendimos dilucidar los genes diana de los elementos que contienen los SNPs. Al analizar los datos de Hi-C de las células HCT116 y centrándonos en los genes expresados diferencialmente (DEG) asociados con el tumor dentro del mismo dominio de asociación topológica (TAD), se identificaron posteriormente algunos genes candidatos, incluidos TAF4, CABLES2, MTG2, RBBP8NL, ADRM1, DIDO1, LINC00659, RPS21, YTHDF1 y HAR1A. Para mapear con precisión su relación regulatoria, generamos dos líneas celulares HCT116 independientes con deleción (KO) utilizando CRISPR-Cas9, cada una con un fragmento corto (89 pb y 136 pb, respectivamente) que abarca el SNP diana rs6061231 eliminado del genoma. Al evaluar el impacto de los genes candidatos anteriores, solo TAF4 y LINC00659 mostraron cambios significativos en la expresión en ambas líneas KO (Fig. S2E). Nuestros datos demostraron que la región alrededor de rs6061231 confiere actividad del potenciador para activar la expresión del gen diana (Fig. 3B), mientras que la deleción del elemento condujo a la regulación ascendente de LINC00659, excluyendo la relación regulatoria directa entre ellos. En contraste, la expresión de TAF4 se redujo notablemente (P < 0.001) tras la deleción del potenciador (Fig. S2F), lo que apoya la posibilidad de que TAF4 fuera el gen diana real del SNP rs6061231.

Para consolidar aún más la regulación específica del alelo de TAF4 por el SNP rs6061231, generamos genotipos HCT116 homocigotos C/C y A/A mediante la edición de CRISPR-Cas9. El análisis de qPCR mostró una expresión significativamente elevada de TAF4 (P < 0.001) en las células homocigotas C/C en comparación con las células A/A. De manera similar, las células homocigotas rs67941642 T/T también exhibieron una expresión significativamente mayor de TAF4 que las células C/C (P < 0.05) (Fig. 3F). Estos resultados fueron muy consistentes con el ensayo de actividad del potenciador (Fig. 3B), lo que demostró la regulación específica del alelo de TAF4 por rs6061231 y rs67941642. La expresión de LINC00659, MTG2, DIDO1, RBBP8NL y OSBPL2 no se alteró significativamente por los cambios de un solo nucleótido (evaluado mediante una prueba t pareada, Fig. S2G). Este resultado, combinado con nuestros datos de secuenciación ChIP-seq de H3K27ac que muestran una mayor actividad del potenciador en el locus LINC00659 tras la deleción (Fig. S2H), indica que el potenciador regula directa y específicamente TAF4 y que la regulación ascendente de LINC00659 es una consecuencia indirecta de las perturbaciones estructurales en lugar de una relación regulatoria directa. Para evaluar el papel de los factores de transcripción GFI1 y p53 unidos a los dos SNPs en la regulación de la expresión de TAF4, realizamos experimentos de desactivación dirigidos a estos factores de transcripción. La depleción de GFI1 o p53 redujo significativamente la expresión de TAF4 (Fig. 3G). Estos resultados establecieron que los potenciadores que albergan estos SNPs activaron la transcripción de TAF4 mediante el reclutamiento de GFI1 y p53.

Estudios previos han encontrado que algunos genes exhiben un papel epistático sobre otros genes en la expresión del fenotipo. Aquí, revelamos elementos reguladores epistáticos que mostraron un papel dominante sobre otros potenciadores en la regulación de la transcripción del gen diana. El resultado subrayó aún más la importancia de analizar exhaustivamente la función molecular de las variantes en desequilibrio de enlace y adyacentes a los resultados de GWAS.

Estudios previos han demostrado que numerosas variantes genéticas no codificantes funcionales pueden alterar la unión de los factores de transcripción. El análisis de motivos sugirió la unión específica de los factores de transcripción GFI1 a rs67941642-T, p53 a rs6061231-C y SIX homeobox 2 (SIX2) a rs6061231-A (Fig. 3D y Fig. S2A). La genotipificación confirmó la heterocigosidad de ambos SNPs rs67941642 y rs6061231 en las células HCT116 y SW480 (Fig. S2B). Posteriormente, la inmunoprecipitación de cromatina (ChIP) seguida de la reacción en cadena de la polimerasa cuantitativa (qPCR) y la secuenciación del amplicón ChIP reveló el reclutamiento preferencial de GFI1 a rs67941642-T (P < 0.001) y de p53 a rs6061231-C (P < 0.001), mientras que la unión de SIX2 no mostró un sesgo alélico significativo en rs6061231 (Fig. 3E y Fig. S2C). Dada los efectos contrastantes de ambos SNPs en la actividad del potenciador, investigamos si los dos factores de transcripción podrían tener un papel interferente en este locus. Por lo tanto, desactivamos individualmente p53 y GFI1 en las células HCT116 y encontramos que la unión de GFI1 en rs67941642 o de p53 en rs6061231 no se vio significativamente afectada, respectivamente, lo que sugiere que no existe una influencia mutua entre estos dos factores de transcripción (Fig. S2D).

SNP-STARR-seq descubrió SNPs que potencialmente contribuyen a la metástasis del CRC

La metástasis es una de las principales causas de muerte en pacientes con CRC, impulsada por complejas alteraciones moleculares adquiridas durante la progresión tumoral. Para identificar sistemáticamente variantes genéticas funcionales que contribuyen a la transición metastásica, aprovechamos un modelo celular isogénico único, comparando la actividad regulatoria de los SNPs en SW480 (una línea celular epitelial derivada de un CRC primario) y SW620 (derivada de un tumor metastásico del ganglio linfático del mismo paciente con CRC que SW480). Este modelo celular pareado proporcionó un sistema excepcionalmente potente para identificar los cambios genéticos y epigenéticos adquiridos en la metástasis, minimizando las diferencias genéticas interindividuales que podrían confundir los resultados. Hipotetizamos que los loci que exhiben función en la modulación de la actividad del potenciador específicamente en las células metastásicas SW620, pero no en las células SW480, representan posibles factores genéticos o epigenéticos asociados con los rasgos metastásicos.

Utilizando el sistema SNP-STARR-seq, analizamos la presencia de SNPs que exhiben diferencias alelo-específicas en la actividad de los potenciadores entre las líneas celulares SW480 y SW620 (tabla S3). Los paSNPs identificados en las células SW620 fueron predominantemente variantes funcionales, con un 88,9% respaldado por evidencia de eQTL y un 72,3% que se predijo que residían dentro de sitios de unión a factores de transcripción (TF) (fig. S3A). La comparación cruzada de los resultados de las dos líneas celulares identificó 3136 SNPs que exhibieron actividad de potenciador alelo-específica específicamente en las células SW620 (fig. 4A y tabla S6). El análisis HOMER predijo un enriquecimiento significativo de motivos de unión para los TFs de la familia bZIP alrededor de estos SNPs (fig. 4, B y C), una familia estructural de TFs predominantemente asociada con la metástasis del CCRC (^[16], ^[22]–^[25]).

Luego, buscamos refinar aún más el mapeo de variantes genéticas funcionales y loci que contribuyen a la metástasis del CCRC. El análisis de RNA-seq que comparó las células SW480 y SW620 identificó 2560 genes que se sobreexpresaron de manera significativa y específica en SW620, lo que sugiere su fuerte asociación con la metástasis del CCRC [P ajustada < 0,05, cambio de pliegue logarítmico (log2FC) > 0,75; fig. S3B y tabla S7]. Los análisis de enriquecimiento funcional [Ontología de genes (GO) y Enciclopedia de genes y genomas de Kioto (KEGG)] confirmaron que estos genes sobreexpresados estaban estrechamente relacionados con procesos relacionados con la metástasis, como la migración celular y las vías de señalización de Wnt (fig. S3C). Entre los 3136 SNPs reguladores específicos de SW620, 2517 (80,3%) se ubicaron dentro de ±1 Mbp de los sitios de inicio de la transcripción de estos genes sobreexpresados en relación con la metástasis, lo que sugiere una relación regulatoria entre las variantes genéticas y estos genes. Al integrar estos hallazgos con los picos de ChIP-seq de H3K27ac específicos de SW620, demostramos que 1697 de estos SNPs se situaron en potenciadores activados por la metástasis. En conjunto, estos resultados subrayan que SNP-STARR-seq es una herramienta eficaz para identificar variantes genéticas funcionales de manera específica para cada tipo de célula, lo que facilita el descubrimiento de SNPs de riesgo asociados con la metástasis del CCRC.

El papel del SNP rs1962004 en la regulación de la expresión de LRRC61

Dado el marcado enriquecimiento de los motivos de unión de los TFs de la familia bZIP en estos loci, hipotetizamos que los SNPs que interrumpen la unión de bZIP impulsan el secuestro de potenciadores durante la progresión metastásica. Por lo tanto, priorizamos 22 SNPs asociados con la metástasis de alta confianza (fig. 4D y tabla S6) con base en cuatro criterios, que incluyen: (i) actividad de potenciador alelo-específica significativa en SW620 (pero negativa en SW480), (ii) ubicación dentro de los picos de H3K27ac específicos de SW620, (iii) proximidad (±1 Mbp) a los genes sobreexpresados en relación con la metástasis y (iv) capacidad para interrumpir los supuestos sitios de unión de los TFs de la familia bZIP. Entre estas variantes priorizadas, encontramos que un SNP, rs1962004, presentaba diferentes genotipos entre las células SW480 y SW620 (SW480, G/G; y SW620, G/A; fig. 4E), lo que indica que el alelo "A" específico de SW620 era muy probablemente una mutación adquirida somáticamente y podría contribuir a la metástasis de las células del CCRC. Por lo tanto, seleccionamos el SNP para una posterior caracterización mecanística. En comparación con las células SW480, las células SW620 mostraron un enriquecimiento adquirido de H3K27ac en este locus, lo que revela la actividad del potenciador de la vecindad. Mientras tanto, observamos una expresión significativamente más fuerte del gen adyacente LRRC61 en las células SW620 en comparación con las células SW480, lo que indica la posible relación regulatoria entre el potenciador específico de la metástasis y la expresión del gen LRRC61 (fig. 4F).

Como se discutió anteriormente, las mutaciones no codificantes generalmente conducen a una unión alterada de los TFs y a una expresión aberrante del gen diana. Se predijo que el SNP rs1962004 residía dentro de un sitio de unión putativo de los TFs de la familia bZIP, JUN y JUND. Por lo tanto, primero confirmamos la unión significativa de JUN y JUND a este locus utilizando ChIP-qPCR (fig. 4G). La redundancia funcional se alinea con los mecanismos establecidos de los TFs bZIP, donde JUN y JUND forman con frecuencia heterodímeros como parte del complejo AP-1 para regular cooperativamente los genes diana involucrados en la oncogénesis y la metástasis (^[26]–^[28]). En particular, estos TFs prefirieron el alelo "A" adquirido por SW620, lo que sugiere que se creó un nuevo sitio de unión de JUN y JUND en el locus durante la metástasis cuando SW620 experimentó una transición de "G" a "A" de rs1962004. Esto también fue consistente con los cambios epigenéticos observados en el locus circundante de las células metastásicas de SW480 a SW620 (fig. 4F). Además, un ensayo de reportero de luciferasa demostró una actividad significativamente mayor del fragmento genómico que contiene el alelo "A" en comparación con el alelo "G" (P < 0,001, fig. 4H).

Además, utilizando la edición del genoma basada en CRISPR-Cas9, generamos dos clones celulares con genotipos SNP rs1962004 homocigotos: G/G y A/A, respectivamente. Nuestros resultados mostraron de manera consistente que las células con el genotipo A/A exhibieron una capacidad proliferativa y migratoria significativamente mayor en comparación con las células SW620 con el genotipo rs1962004 heterocigoto (es decir, G/A), que, a su vez, exhibieron una mayor capacidad pro-metastásica que las células homocigotas G/G (fig. 4, I a K). Complementariamente, la deleción de la región del potenciador de 104 pares de bases que abarca rs1962004 (fig. S3, D y E) atenuó significativamente la proliferación en las células SW620. Mientras tanto, la expresión de LRRC61 se redujo significativamente en las células con el potenciador eliminado en comparación con los controles de tipo salvaje (P < 0,01, fig. S3F). En conjunto, estos hallazgos demuestran que tanto el alelo "A" adquirido como el elemento del potenciador en sí impulsan mecánicamente la agresividad metastásica. La unión de JUN o JUND, causada por la transición de "G" a "A" en rs1962004, mejoró posteriormente la actividad cis-reguladora y aumentó la expresión de LRRC61.

El papel metastásico de LRRC61

Para establecer aún más la relación regulatoria entre el SNP y LRRC61, consultamos los datos de eQTL del grupo TCGA-CRC y revelamos que LRRC61 exhibió una expresión correlacionada con los genotipos de rs1962004 (P = 7,19 × 10^19). Es decir, los individuos que portaban el genotipo A/A tenían una expresión significativamente mayor de LRRC61, seguido por los portadores de A/G y G/G (fig. 5A). En consonancia con las tendencias fenotípicas más amplias observadas a nivel de la población, las células SW620 editadas con un genotipo A/A también mostraron niveles elevados de LRRC61 en comparación con las células de la misma procedencia que portaban el genotipo G/G (fig. 5B). El análisis previo reveló un enriquecimiento relacionado con la metástasis de JUN y JUND en el locus rs1962004. Aquí, realizamos la supresión mediada por siRNA de JUN y JUND, ambos de los cuales atenuaron la expresión de LRRC61 (fig. 5, C y D).

La superfamilia de proteínas que contienen repeticiones ricas en leucina (LRRC) tiene funciones complejas y diversas en varios tumores sólidos (como el sarcoma, el melanoma y el glioblastoma), y algunos de sus miembros (como LRRC15 y LRRC32) impulsan la invasión y la metástasis del tumor al regular la señalización del factor de crecimiento transformante-β, promover la activación de fibroblastos asociados con el cáncer e influir en el microambiente inmunitario (^[29], ^[30]). Sin embargo, el papel de LRRC61 en el cáncer, especialmente en el CCRC, no se había comprendido bien. Aquí, pretendimos investigar su posible participación en el impulso de los fenotipos metastásicos. La supresión de LRRC61 (P < 0,001) utilizando ARN de cadena corta (shRNA) condujo a una disminución de la capacidad de formación de colonias y la migración (fig. 5, E a G), que son fenotipos estrechamente asociados con la metástasis. La metástasis causa significativamente la mortalidad del paciente y un mal pronóstico. El seguimiento clínico de los pacientes con CCRC reveló que la expresión de LRRC61 era mayor en los tejidos tumorales en comparación con los tejidos normales (fig. 5H), y los pacientes con una expresión elevada de LRRC61 generalmente tenían un peor pronóstico (fig. 5I). Un modelo genético para esta regulación propone que el alelo de riesgo rs1962004-A, al mejorar la unión de JUN/JUND y aumentar la expresión de LRRC61, subyace al aumento del potencial metastásico y al peor pronóstico del paciente (fig. 5J).

Methyl-STARR-seq identifica elementos potenciadores dependientes de la metilación

Además de las variantes genéticas que influyen en los fenotipos de la enfermedad, las variaciones epigenéticas también se reconocen ampliamente como factores etiológicos clave que afectan la salud humana. La metilación del ADN representa un mecanismo epigenético crítico que regula la expresión génica. Si bien la investigación existente se centra predominantemente en la metilación del promotor, los papeles regulatorios y el mapeo genómico preciso de la metilación dentro de las regiones de los potenciadores, que constituyen una proporción sustancial del epigenoma, siguen estando insuficientemente caracterizados. Para abordar esta brecha, diseñamos una biblioteca de oligonucleótidos dirigida a las regiones de potenciador metiladas diferencialmente en el CCRC y a ventanas de 120 nt centradas en los picos de H3K27ac. Esta biblioteca se sintetizó y clonó en el vector CpG-free Methyl-STARR-seq (^[14]). Tras el tratamiento enzimático in vitro con glicerol de control o M.SssI, una enzima que puede catalizar la metilación completa de la secuencia CpG únicamente en los fragmentos de inserción (^[31]), generamos las versiones no metiladas y completamente metiladas de los mismos grupos de plásmidos. La tasa de conversión y el estado de metilación se confirmaron mediante el secuenciamiento de bisulfito (BS-seq; fig. S4, A y B). Estos constructos se transfectaron en las líneas celulares HCT116, SW480 y SW620, respectivamente, para cuantificar su actividad de potenciador. Posteriormente, se construyeron bibliotecas de ARNm y ADN de entrada y se sometieron a secuenciación y análisis bioinformáticos (esquema del flujo de trabajo experimental en la fig. 1). Al calcular la actividad diferencial del potenciador (EAS) entre los estados metilados y no metilados, identificamos los elementos reguladores más influenciados por la metilación del ADN en cada línea celular (HCT116, n = 1763; SW480, n = 3588; y SW620, n = 4112; fig. S4, C y D, y tabla S4).

Tras la metilación, aproximadamente el 60% de los fragmentos probados exhibieron una actividad de potenciador reducida en las células HCT116, en comparación con aproximadamente el 40% en las células SW480 y SW620, lo que destaca la heterogeneidad específica de la línea celular (fig. 6A). Para los elementos comunes a HCT116 y SW480, la validación utilizando un vector de reportero de Lucia CpG-free (InvivoGen) demostró una alta concordancia con los resultados de secuenciación de alto rendimiento (fig. 6B y fig. S4E). Las alteraciones dependientes de la metilación en la actividad del potenciador observadas probablemente se debieron a la sensibilidad a la metilación de los TFs. El análisis del potencial de unión de los TFs reveló sitios de unión predichos dentro de más del 95% de estos elementos (tabla S8). En particular, casi el 80% de estos sitios contenían motivos ricos en CpG dentro de sus secuencias de unión de TFs predichas, lo que sugiere que la metilación de CpG podría afectar directamente la ocupación de los TFs (fig. 6C). Además, los TFs enriquecidos en estos elementos reguladores identificados por Methyl-STARR-seq mostraron una fuerte concordancia con las categorías de TFs methylplus (definidas como la metilación de CpG dentro de los sitios de unión aumenta la afinidad de unión) y methylminus (definidas como la metilación dentro de los sitios de unión de CpG disminuye la afinidad de unión) definidas por Methylation Systematic Evolution of Ligands by Exponential Enrichment (Methyl-SELEX) (^[32]), lo que confirma la fiabilidad de nuestro sistema de detección y también el papel de los TFs en la modulación de la actividad del potenciador dependiente de la metilación (fig. 6D). Realizamos además un análisis de enriquecimiento de motivos HOMER exhaustivo para los elementos methylminus frente a los elementos methylplus identificados en las líneas celulares HCT116, SW480 y SW620, así como específicamente para los 3008 elementos reguladores sensibles a la metilación adquiridos por la metástasis únicos de SW620. Nuestros resultados revelaron que los elementos methylminus en las líneas celulares HCT116 y SW620 más agresivas se enriquecieron preferentemente con los TFs de la familia ETS, que son impulsores bien establecidos de la invasión y la metástasis (^[33], ^[34]). En contraste, la línea primaria SW480 mostró un mayor enriquecimiento de los motivos bZIP/AP-1. Este notable patrón específico del contexto, en particular el fuerte enriquecimiento del factor ETS en los elementos sensibles a la metilación adquiridos por la metástasis de SW620, proporciona un vínculo mecanístico convincente entre la reprogramación epigenética durante la progresión metastásica y la activación de los programas transcripcionales pro-metastásicos (fig. S5A).

Posteriormente, nos centramos en las DMR (regiones de metilación diferencial) que podrían desempeñar un papel fundamental en la patogénesis del CCRC. Integrando el análisis de motivos TF (factores de transcripción) con los datos de modificación de H3K4me1, identificamos seis sitios CpG funcionalmente putativos, entre ellos cg25982657, cg08640619, cg23516310, cg14302214, cg14179401 y cg14106234, ninguno de los cuales había sido caracterizado previamente (fig. 6E y tabla S9). El análisis de los datos de metilación de muestras tumorales y normales emparejadas de 14 tipos de cáncer reveló que cg25982657 y cg08640619 estaban hipermetilados específicamente en el CCRC (| | > 0,25, tasa de descubrimiento falso (FDR) < 0,05; fig. S5B y tabla S10].

Utilizando cg08640619 como un sitio representativo, realizamos una validación funcional. Este sitio CpG se encontraba en el primer intrón del gen KIRREL1, que exhibía una actividad potenciadora significativa (fig. 6F). La predicción in silico indicó una posible unión de RUNX2, un TF conocido por unirse preferentemente al ADN no metilado. Para confirmar el papel de la metilación del ADN en la unión de RUNX2 a este locus, utilizamos la metilación dirigida mediada por CRISPR-Cas9 inactiva (dCas9) en las células HCT116 y SW480. El ChIP-qPCR posterior demostró una reducción significativa del enriquecimiento de RUNX2 en el locus hipermetilado en las células HCT116 (P < 0,001, fig. 6G). La inhibición de la expresión de RUNX2 provocó una regulación a la baja significativa de sus genes adyacentes, KIRREL1 y ETV3, dentro del mismo TAD (dominio topológicamente asociado), lo que confirmó que RUNX2 se unía preferentemente a la secuencia no metilada y regulaba positivamente KIRREL1 y ETV3 (fig. 6H). Tras observar una correlación negativa significativa entre la metilación de cg08640619 y la expresión de KIRREL1/ETV3 en muestras clínicas (fig. 6I), buscamos validar funcionalmente este hallazgo. Utilizando la hipermetilación dirigida mediada por CRISPR-dCas9-DNMT3A, confirmamos que el aumento de la metilación causa directamente una regulación a la baja significativa de ambos genes, lo que demuestra una relación regulatoria causal (fig. S5C). Se encontró que la expresión de KIRREL1 estaba significativamente aumentada en los tejidos tumorales en comparación con los tejidos normales correspondientes en las muestras emparejadas de adenocarcinoma de colon (COAD) del TCGA (fig. S5D). Este aumento de la expresión puede estar asociado con el estado de hipometilación de cg08640619 en los tumores de CCRC. El análisis pronóstico de los datos clínicos reveló que la hipometilación de cg08640619 está significativamente asociada con un mal pronóstico del paciente en la enfermedad en estadio avanzado (estadios III a IV) (fig. 6J).

Dado los patrones de hipometilación específicos del cáncer (fig. S5B), hipotetizamos que cg08640619 y cg25982657 podrían servir como biomarcadores sensibles para la detección temprana del CCRC. Para evaluar el potencial diagnóstico de los sitios CpG identificados, aplicamos tanto la regresión logística (LogR) como los clasificadores de máquinas de vectores de soporte (SVM) a muestras de 14 tipos de cáncer diferentes. Además, se realizó una validación cruzada de 10 pliegues para garantizar la robustez del modelo (fig. S5, E, G y H). Los resultados mostraron que ambos sitios CpG discriminaban eficazmente los tumores de CCRC de los tejidos normales, pero no tan bien para otros tumores (fig. 6K y fig. S5F). Los variantes epigenéticos funcionales detectados por Methyl-STARR-seq estaban altamente correlacionados con el riesgo de enfermedad, lo que les confiere un gran potencial clínico para el diagnóstico específico y sensible del cáncer.

De manera similar a nuestro cribado de SNPs funcionales asociados con la metástasis, identificamos 3008 regiones reguladoras sensibles a la metilación que se desregulaban específicamente en las células SW620 (fig. 6L y tabla S11). El análisis integrativo de los perfiles de modificación de H3K27ac y los datos de RNA-seq en ambas líneas celulares identificó 182 fragmentos genómicos que exhibían una actividad reguladora mejorada en el estado de hipometilación (fig. 6M). Estos loci mostraron simultáneamente un enriquecimiento significativamente elevado de H3K27ac en las células SW620 en comparación con las células SW480 y se ubicaron dentro de ± 1 Mb de genes significativamente sobreexpresados en esta línea celular metastásica.

Tomemos como ejemplo un sitio putativo, cg02205193. Estos dinucleótidos CpG se encuentran cerca del promotor de LINC00460, un gen que codifica un ARN endógeno competidor que, según se ha informado en varios estudios, promueve la metástasis del CCRC (^[35]–^[37]). La región genómica que encierra cg02205193 exhibió señales de H3K27ac significativamente más altas en las células SW620 en comparación con las células SW480 (fig. 6N). En consecuencia, LINC00460 está significativamente sobreexpresado en los pacientes con CCRC metastásico, como se observó en la línea celular metastásica SW620, lo que podría atribuirse a su estado de hipometilación (^[38]). Además, la metilación en cg02205193 influyó supuestamente en la unión de GFI1, un TF (factor de transcripción) que no se une a ADN metilado, según los datos de Methyl-SELEX (^[32]). Específicamente, el estado de hipometilación de cg02205193 facilita una mayor captación de GFI1, lo que, a su vez, activa la expresión del gen diana LINC00460 y promueve la metástasis del CCRC. Basándonos en Hi-C y la expresión diferencial de genes entre las células SW480/SW620, también identificamos regiones reguladoras que encierran sitios CpG previamente no caracterizados cuyo estado de metilación se asoció con la expresión de los genes TCF12 y HNF1, dos genes establecidos implicados en la metástasis del CCRC (fig. S5I) (^[39], ^[40]).

Por lo tanto, al integrar el cribado funcional de Methyl-STARR-seq con la validación multiómica, pudimos identificar eficazmente posibles biomarcadores diagnósticos. Este enfoque no solo delimita un paisaje regulatorio preciso de la disregulación epigenética en el CCRC, sino que también proporciona una base mecánica para la identificación de biomarcadores diagnósticos y terapéuticos previamente no identificados con una amplia aplicabilidad clínica.

Diseño del proyecto para identificar variantes funcionales en el CCRC primario y metastásico

Cómo las variantes reguladoras dependientes del contexto impulsan la patogénesis del CCRC, particularmente durante la metástasis, no se había estudiado adecuadamente en estudios anteriores. Aquí, desarrollamos una estrategia que pudo evaluar sistemáticamente tanto las variantes genéticas como las epigenéticas para su impacto en la actividad reguladora transcripcional de los elementos potenciadores del CCRC (fig. 1).

En primer lugar, definimos un conjunto de posibles regiones potenciadoras del CCRC identificando 18.880 intervalos genómicos que se superponían entre los picos de ChIP-seq de H3K27ac y los picos de ATAC-seq en la línea celular de CCRC HCT116, lo que indica tanto una actividad reguladora activa como la accesibilidad de la cromatina. Para distinguir el paisaje regulatorio del CCRC de los elementos pancancerosos comunes, realizamos un análisis epigenómico comparativo en 20 tipos de cáncer (tablas S1 y S2), lo que nos permitió clasificar estos potenciadores como comunes en muchos cánceres o más restringidos en su perfil de actividad.

Posteriormente, utilizamos dos sistemas distintos, SNP-STARR-seq y Methyl-STARR-seq, para evaluar la actividad del potenciador que contenía las variantes genéticas y epigenéticas, respectivamente. Para la biblioteca SNP-STARR-seq, sintetizamos una biblioteca genómica de 120 nucleótidos (nt) que abarcaba todos los sitios SNP dentro de las 18.880 regiones potenciadoras putativas y los loci SNP de éxito en el GWAS (estudio de asociación del genoma completo) y los loci SNP vinculados al GWAS. Para cada SNP, se incluyó una secuencia genómica de 120 pb centrada en la variante. Se sintetizaron todas las posibles variantes alélicas (A, T, C y G) en el sitio polimórfico como oligonucleótidos individuales, lo que permitió una evaluación exhaustiva de los efectos alélicos específicos de las variantes genéticas humanas conocidas. Para la biblioteca de metilación, nos centramos en dos categorías de secuencias genómicas: (i) ventanas de 120 nt centradas en los picos de H3K27ac, ya que estas regiones potenciadoras y promotoras activas son puntos calientes establecidos para la reprogramación epigenética en el cáncer, donde los cambios en la metilación del ADN ejercen su impacto funcional más profundo. (ii) regiones genómicas que corresponden a las DMR del CCRC. Tras la construcción de la biblioteca de plásmidos, realizamos un tratamiento enzimático in vitro utilizando la metiltransferasa M.SssI para generar dos bibliotecas distintas: biblioteca de plásmidos totalmente metilada (tratada con M.SssI) y biblioteca de plásmidos de control no metilada (tratada con un simulacro) (ver Materiales y Métodos). Las puntuaciones de actividad del potenciador (EAS) se derivaron de los recuentos de ARN/ADN normalizados de DESeq2 como EAS = (ARN + 1)/ADN. Los cambios de alelo o metilación se calcularon comparando los valores de EAS entre el alelo alternativo y el alelo de referencia o la metilación y la no metilación, con significación estadística. Mediante un cribado sistemático en tres líneas celulares de CCRC, modelos de adenocarcinoma primarios (HCT116 y SW480) y una línea celular derivada metastásica (SW620), identificamos variantes reguladoras funcionales moduladas significativamente tanto por las mutaciones genéticas como por las alteraciones epigenéticas. Los SNPs preferentemente activos (paSNPs; definidos como SNPs que confieren una puntuación de potenciador diferencial significativa y un P ajustado < 0,05; ver Materiales y Métodos) y los sitios de metilación preferentemente activos (paMethylsites; definidos como la metilación de CpG que confiere una puntuación de potenciador diferencial significativa y un P ajustado < 0,05; ver Materiales y Métodos) identificados en las tres líneas celulares se presentan en las tablas S3 y S4.

Al integrar aún más los conjuntos de datos clínicos de RNA-seq, eQTL y epigenómicos, asignamos las variantes reguladoras a sus genes diana. Centrándonos en las variantes implicadas en el desarrollo y la metástasis del CCRC, validamos la relevancia funcional de las variantes identificadas a través de esta vía de cribado utilizando técnicas que incluyen ensayos de reportero de genes, inmunoprecipitación de cromatina (ChIP) y edición genómica basada en CRISPR. En conjunto, este marco de genómica funcional integrado no solo descifra los vínculos mecánicos entre las variantes no codificantes y la patogénesis/metástasis del CCRC, sino que también identifica previamente biomarcadores diagnósticos y terapéuticos no identificados con un potencial de traslación significativo.

Dentro de los 18.880 potenciadores candidatos, seleccionamos los SNPs con una frecuencia alélica menor (MAF) ≥ 0,05 tanto en las poblaciones europeas como en las asiáticas orientales del Proyecto de los 1000 Genomas (^[15]), lo que resultó en 24.526 SNPs comunes. Para investigar el potencial regulador de los SNPs dentro de las regiones del potenciador, incorporamos datos del Catálogo de GWAS y recuperamos 238 loci asociados al CCRC significativos a nivel del genoma, junto con 6125 SNPs adicionales en LD [coeficiente de determinación (R2) ≥ 0,5] con estos loci. Al combinar los SNPs identificados a partir de las regiones del potenciador del CCRC con los SNPs relacionados con el GWAS, obtuvimos 30.790 SNPs únicos como elementos reguladores candidatos para el análisis de cribado funcional. Para cada SNP, se extrajo una secuencia genómica de 120 pb centrada en la variante (± 60 pb). Se sintetizaron todas las variantes alélicas como oligonucleótidos individuales. En total, diseñamos y sintetizamos 91.566 oligonucleótidos únicos, formando una biblioteca de SNP-STARR-seq de CCRC para un ensayo de reportero paralelo masivo (fig. 2A). La agrupación jerárquica demostró una alta reproducibilidad en ambos replicados biológicos y técnicos [reacción en cadena de la polimerasa (PCR)] del mismo tipo de célula (fig. 2B). Para cada SNP, el efecto alélico específico en la actividad del potenciador se cuantificó calculando el cambio transformado logarítmicamente de EAS entre el alelo alternativo y el alelo de referencia. Todos los SNPs calificados se clasificaron en función de este valor y se seleccionaron los SNPs con actividades de potenciador diferencialmente significativas para la posterior caracterización funcional (HCT116, n = 1969 SNPs; y SW480, n = 1846 SNPs; fig. 2C y tabla S3). La validación de los resultados de SNP-STARR-seq utilizando un ensayo de reportero de luciferasa de un solo locus demostró una fuerte concordancia con los resultados de secuenciación de alto rendimiento en ambas líneas celulares (fig. 2D y fig. S1A).

Tanto las líneas celulares HCT116 como SW480 representan el adenocarcinoma colorrectal en etapa temprana, que exhibió una superposición significativa de variantes funcionales (P = 0.048, prueba exacta de Fisher; Fig. 2E). A continuación, nos centramos en 922 SNPs comunes que exhiben efectos significativos en ambas líneas celulares (Tabla S5). La función de estas variantes no codificantes también está respaldada por datos de eQTL (82.65%), lo que demuestra su impacto regulador en la expresión génica proximal o distal (Fig. 2F). Además, se predijo computacionalmente que más del 75% de estas variantes albergan sitios de unión a factores de transcripción (Fig. 2G). Se obtuvieron resultados similares de los análisis de cada línea celular individualmente (Fig. S1, B y C). El análisis de motivos mediante la optimización hipergeométrica de enriquecimiento de motivos (HOMER) reveló además un enriquecimiento significativo de factores de transcripción (FT) asociados con el CRC dentro de los segmentos genómicos que muestran diferencias significativas en la actividad del potenciador alélico, incluidos el factor de transcripción tipo factor de cremallera leucina básica (BATF), FOS like 2, la subunidad del factor de transcripción AP-1 (FOSL2) y el antígeno relacionado con Fos-2 (FRA2), que se han relacionado directamente con la progresión y metástasis del CRC (^[16], ^[17]), junto con JUNB, un componente del complejo AP-1 que se sabe que está desregulado en el CRC (Fig. 2H y Fig. S1D) (^[18]). Basándonos en nuestra anotación de potenciadores pancancerosos mencionada anteriormente, el 81.6% de los 922 SNPs comunes residían dentro de las regiones de potenciadores putativos del CRC, aunque algunas variantes vinculadas a estudios de asociación del genoma completo (GWAS) se encuentran fuera de estas regiones. Este enriquecimiento es significativamente mayor que el observado en otros 19 tipos de cáncer, lo que confirma que nuestro sistema de alto rendimiento identificó eficazmente variantes regulatorias relevantes para el contexto del CRC (Fig. 2I). En conjunto, estos hallazgos demostraron que el SNP-STARR-seq podría capturar con éxito la gran mayoría de las variantes funcionales.

El papel epistático del SNP rs67941624 en la modulación de la actividad del potenciador en las células de CRC

Entre los 922 SNPs putativamente funcionales, rs6061231 se encuentra en un locus de riesgo de CRC y se identificó como un SNP principal por los estudios GWAS (^[19]), mientras que el SNP rs67941642 está en completo desequilibrio de enlace (R2 = 1, https://ldlink.nih.gov/) con rs6061231 (Fig. 3A). Ambos SNPs residían dentro de las regiones genómicas que presentaban fuertes marcas epigenéticas asociadas con potenciadores (H3K4me1 y H3K27ac) en las células de CRC humanas (Fig. 3B). Dada la estrecha proximidad y el vínculo genético, utilizamos un ensayo de reportero de luciferasa dual para investigar el papel sinérgico de las dos variantes. Primero, confirmamos la actividad del potenciador de los fragmentos genómicos que contienen cualquiera de los alelos mayor o menor de los SNPs rs6061231 y rs67941642. En particular, el alelo rs67941642-T exhibió una actividad del potenciador significativamente mayor que el alelo C, y el alelo rs6061231-C mostró una actividad significativamente mayor que el alelo A (Fig. 3B), lo que es consistente con el resultado de SNP-STARR-seq (Tabla S5). En la cohorte combinada de Asia Oriental y Europa (EAS + EUR), el análisis de haplotipos del locus rs6061231/rs67941642 reveló frecuencias del 77.0% para el haplotipo C-C (rs6061231-C/rs67941642-C) y del 23.0% para el haplotipo A-T (rs6061231-A/rs67941642-T), mientras que los otros haplotipos (es decir, C-T y A-C) prácticamente no ocurren en la población. Dado que los dos SNPs en desequilibrio de enlace exhibieron actividades de potenciador opuestas, preguntamos qué SNP fue dominante en la determinación de la actividad del potenciador alélico. Por lo tanto, construimos los plásmidos reporteros que contienen los haplotipos C-C y A-T, respectivamente, y los transfectamos a las células de CRC para una comparación cuantitativa. El resultado del ensayo mostró que el haplotipo A-T exhibió una actividad del potenciador significativamente mayor en comparación con el haplotipo C-C (P < 0.001, Fig. 3C), lo que indicó que el SNP rs67941642 confirió un papel predominante en la modulación de la actividad del potenciador dentro de su contexto genómico. El SNP rs6061231 fue un SNP principal en los estudios GWAS, lo que sugiere que el SNP principal resultó ser un proxy de la variante funcional del SNP rs67941642. Si bien ambos SNPs exhibieron actividad del potenciador específica del alelo, el haplotipo A-T (rs6061231-A/rs67941642-T) confirió una actividad significativamente mayor que el haplotipo C-C (Fig. 3C). Esto indica que la producción regulatoria neta de este locus está dominada por el efecto combinado de ambos SNPs, con rs67941642-T desempeñando un papel predominante en la determinación de la fuerza general del potenciador de los haplotipos humanos comunes. Se ha establecido que algunos genes exhiben un papel epistático sobre otros genes en la expresión del fenotipo. Aquí, revelamos elementos reguladores epistáticos que mostraron un papel dominante sobre otros potenciadores en la regulación de la transcripción del gen diana. El resultado subrayó aún más la importancia de analizar exhaustivamente la función molecular de las variantes en desequilibrio de enlace y adyacentes a los resultados de los estudios GWAS.

Estudios previos han demostrado que numerosas variantes genéticas no codificantes funcionales pueden alterar la unión de los factores de transcripción. El análisis de motivos sugirió la unión específica de los factores de transcripción GFI1 a rs67941642-T, p53 a rs6061231-C y SIX homeobox 2 (SIX2) a rs6061231-A (Fig. 3D y Fig. S2A). La genotipificación confirmó la heterocigosidad de ambos SNPs rs67941642 y rs6061231 en las células HCT116 y SW480 (Fig. S2B). Posteriormente, la inmunoprecipitación de cromatina (ChIP) seguida de la reacción en cadena de la polimerasa cuantitativa (qPCR) y la secuenciación del amplicón de ChIP reveló el reclutamiento preferencial de GFI1 a rs67941642-T (P < 0.001) y de p53 a rs6061231-C (P < 0.001), mientras que la unión de SIX2 no mostró un sesgo alélico significativo en rs6061231 (Fig. 3E y Fig. S2C). Dada los efectos contrastantes de ambos SNPs en la actividad del potenciador, investigamos si los dos factores de transcripción podrían tener un papel interferente en este locus. Por lo tanto, desactivamos individualmente p53 y GFI1 en las células HCT116 y encontramos que la unión de GFI1 en rs67941642 o de p53 en rs6061231 no se vio significativamente afectada, respectivamente, lo que sugiere que no existe una influencia mutua entre estos dos factores de transcripción (Fig. S2D).

El locus rs6061231/rs67941642 regula cis la expresión de TAF4

Un estudio previo demostró que la amplificación del número de copias del cromosoma 20q sirve como biomarcador pronóstico en el subtipo de CRC con microsatélites estables (^[20]). Dado que los SNPs rs6061231 y rs67941642 se encuentran dentro del locus Chr20q, nos propusimos comprender su relevancia para la significación clínica. Primero, pretendimos dilucidar los genes diana de los elementos que contienen los SNPs. Al analizar los datos de Hi-C de las células HCT116 y centrándonos en los genes expresados diferencialmente (DEG) asociados con el tumor dentro del mismo dominio de asociación topológica (TAD), se identificaron posteriormente algunos genes candidatos, incluidos TAF4, CABLES2, MTG2, RBBP8NL, ADRM1, DIDO1, LINC00659, RPS21, YTHDF1 y HAR1A. Para mapear con precisión su relación regulatoria, generamos dos líneas celulares HCT116 independientes con deleción (KO) utilizando CRISPR-Cas9, cada una con un fragmento corto (89 pb y 136 pb, respectivamente) que abarca el SNP diana rs6061231 eliminado del genoma. Al evaluar el impacto de los genes candidatos anteriores, solo TAF4 y LINC00659 mostraron cambios significativos en la expresión en ambas líneas KO (Fig. S2E). Nuestros datos demostraron que la región alrededor de rs6061231 confiere actividad del potenciador para activar la expresión del gen diana (Fig. 3B), mientras que la deleción del elemento condujo a la regulación ascendente de LINC00659, excluyendo la relación regulatoria directa entre ellos. En contraste, la expresión de TAF4 se redujo notablemente (P < 0.001) tras la deleción del potenciador (Fig. S2F), lo que respalda la posibilidad de que TAF4 fuera el objetivo real del SNP rs6061231.

Para consolidar aún más la regulación específica del alelo de TAF4 por el SNP rs6061231, generamos genotipos HCT116 homocigotos C/C y A/A mediante la edición de CRISPR-Cas9. El análisis de qPCR mostró una expresión significativamente elevada de TAF4 (P < 0.001) en las células homocigotas C/C en comparación con las células A/A. De manera similar, las células homocigotas rs67941642 T/T también exhibieron una expresión significativamente mayor de TAF4 que las células C/C (P < 0.05) (Fig. 3F). Estos resultados fueron muy consistentes con el ensayo de actividad del potenciador (Fig. 3B), lo que demostró la regulación específica del alelo de TAF4 por rs6061231 y rs67941642. La expresión de LINC00659, MTG2, DIDO1, RBBP8NL y OSBPL2 no se alteró significativamente mediante los cambios de un solo nucleótido (evaluado mediante una prueba t pareada, Fig. S2G). Este resultado, combinado con nuestros datos de secuenciación de ChIP-seq de H3K27ac que muestran una mayor actividad del potenciador en el locus LINC00659 tras la deleción (Fig. S2H), indica que el potenciador regula directa y específicamente TAF4 y que la regulación ascendente de LINC00659 es una consecuencia indirecta de las perturbaciones estructurales en lugar de una relación regulatoria directa. Para evaluar el papel de los factores de transcripción GFI1 y p53 unidos a los dos SNPs en la regulación de la expresión de TAF4, realizamos experimentos de desactivación dirigidos a estos factores de transcripción. La depleción de GFI1 o p53 redujo significativamente la expresión de TAF4 (Fig. 3G). Estos resultados establecieron que los potenciadores que albergan estos SNPs activaron la transcripción de TAF4 mediante el reclutamiento de GFI1 y p53.

Estudios previos han encontrado que numerosos variantes genéticas no codificantes funcionales pueden alterar la unión de los factores de transcripción. El análisis de motivos sugirió la unión específica de los factores de transcripción GFI1 a rs67941642-T, p53 a rs6061231-C y SIX homeobox 2 (SIX2) a rs6061231-A (Fig. 3D y Fig. S2A). La genotipificación confirmó la heterocigosidad de ambos SNPs rs67941642 y rs6061231 en las células HCT116 y SW480 (Fig. S2B). Posteriormente, la inmunoprecipitación de cromatina (ChIP) seguida de la reacción en cadena de la polimerasa cuantitativa (qPCR) y la secuenciación del amplicón de ChIP reveló el reclutamiento preferencial de GFI1 a rs67941642-T (P < 0.001) y de p53 a rs6061231-C (P < 0.001), mientras que la unión de SIX2 no mostró un sesgo alélico significativo en rs6061231 (Fig. 3E y Fig. S2C). Dada los efectos contrastantes de ambos SNPs en la actividad del potenciador, investigamos si los dos factores de transcripción podrían tener un papel interferente en este locus. Por lo tanto, desactivamos individualmente p53 y GFI1 en las células HCT116 y encontramos que la unión de GFI1 en rs67941642 o de p53 en rs6061231 no se vio significativamente afectada, respectivamente, lo que sugiere que no existe una influencia mutua entre estos dos factores de transcripción (Fig. S2D).

El locus rs6061231/rs67941642 regula cis la expresión de TAF4

Un estudio previo demostró que la amplificación del número de copias del cromosoma 20q sirve como biomarcador pronóstico en el subtipo de CRC con microsatélites estables (^[20]). Dado que los SNPs rs6061231 y rs67941642 se encuentran dentro del locus Chr20q, nos propusimos comprender su relevancia para la significación clínica. Primero, pretendimos dilucidar los genes diana de los elementos que contienen los SNPs. Al analizar los datos de Hi-C de las células HCT116 y centrándonos en los genes expresados diferencialmente (DEG) asociados con el tumor dentro del mismo dominio de asociación topológica (TAD), se identificaron posteriormente algunos genes candidatos, incluidos TAF4, CABLES2, MTG2, RBBP8NL, ADRM1, DIDO1, LINC00659, RPS21, YTHDF1 y HAR1A. Para mapear con precisión su relación regulatoria, generamos dos líneas celulares HCT116 independientes con deleción (KO) utilizando CRISPR-Cas9, cada una con un fragmento corto (89 pb y 136 pb, respectivamente) que abarca el SNP diana rs6061231 eliminado del genoma. Al evaluar el impacto de los genes candidatos anteriores, solo TAF4 y LINC00659 mostraron cambios significativos en la expresión en ambas líneas KO (Fig. S2E). Nuestros datos demostraron que la región alrededor de rs6061231 confiere actividad del potenciador para activar la expresión del gen diana (Fig. 3B), mientras que la deleción del elemento condujo a la regulación ascendente de LINC00659, excluyendo la relación regulatoria directa entre ellos. En contraste, la expresión de TAF4 se redujo notablemente (P < 0.001) tras la deleción del potenciador (Fig. S2F), lo que respalda la posibilidad de que TAF4 fuera el objetivo real del SNP rs6061231.

Para consolidar aún más la regulación específica del alelo de TAF4 por el SNP rs6061231, generamos genotipos HCT116 homocigotos C/C y A/A mediante la edición de CRISPR-Cas9. El análisis de qPCR mostró una expresión significativamente elevada de TAF4 (P < 0.001) en las células homocigotas C/C en comparación con las células A/A. De manera similar, las células homocigotas rs67941642 T/T también exhibieron una expresión significativamente mayor de TAF4 que las células C/C (P < 0.05) (Fig. 3F). Estos resultados fueron muy consistentes con el ensayo de actividad del potenciador (Fig. 3B), lo que demostró la regulación específica del alelo de TAF4 por rs6061231 y rs67941642. La expresión de LINC00659, MTG2, DIDO1, RBBP8NL y OSBPL2 no se alteró significativamente mediante los cambios de un solo nucleótido (evaluado mediante una prueba t pareada, Fig. S2G). Este resultado, combinado con nuestros datos de secuenciación de ChIP-seq de H3K27ac que muestran una mayor actividad del potenciador en el locus LINC00659 tras la deleción (Fig. S2H), indica que el potenciador regula directa y específicamente TAF4 y que la regulación ascendente de LINC00659 es una consecuencia indirecta de las perturbaciones estructurales en lugar de una relación regulatoria directa. Para evaluar el papel de los factores de transcripción GFI1 y p53 unidos a los dos SNPs en la regulación de la expresión de TAF4, realizamos experimentos de desactivación dirigidos a estos factores de transcripción. La depleción de GFI1 o p53 redujo significativamente la expresión de TAF4 (Fig. 3G). Estos resultados establecieron que los potenciadores que albergan estos SNPs activaron la transcripción de TAF4 mediante el reclutamiento de GFI1 y p53.

Estudios previos han encontrado que el Taf4 murino antagoniza la inhibición mediada por PRC2 de la expresión de los genes de las células madre para mantener la homeostasis intestinal (^[21]). Para sondear la función supresora de tumores de TAF4, realizamos la desactivación de TAF4 en las células HCT116 utilizando ARN de interferencia pequeño (siRNA; Fig. S2I). La falta de TAF4 promovió significativamente tanto la proliferación celular (Fig. 3H) como la migración (Fig. 3I), lo que refleja los fenotipos malignos. Para relacionar mecánicamente esta función de supresión del crecimiento directamente con la actividad del potenciador, analizamos la proliferación de las células KO del potenciador. Coherente con la pérdida de un regulador del crecimiento inhibidor, estas células KO del potenciador exhibieron una mayor proliferación celular (Fig. S2J). Además, evaluamos las consecuencias fenotípicas de los efectos reguladores específicos del alelo realizando ensayos de proliferación en las líneas celulares homocigotas con mutación puntual. Para rs67941642, las células con el genotipo T/T (asociado con una mayor expresión de TAF4) mostraron una tendencia a una menor proliferación en comparación con las células C/C, lo que se alinea con el papel supresor de tumores de TAF4. En contraste, para rs6061231, el cambio de genotipo de A/A a C/C no produjo una diferencia estadísticamente significativa en el crecimiento celular.

El análisis de supervivencia utilizando los datos del Programa del Genoma del Cáncer (TCGA) reveló que los pacientes con baja expresión (37%) de TAF4 sufrieron una tasa de supervivencia general a 5 años significativamente peor (57%) en comparación con aquellos con alta expresión de TAF4 (P < 0.05, Fig. 3J). Estos hallazgos establecieron que TAF4 es un supresor de tumores bajo la regulación de dos SNPs de riesgo de CRC altamente conectados, rs6061231/rs67941642, y que el alelo con una actividad del potenciador general sinérgicamente menor se asoció con un mal resultado clínico. Como hemos comentado anteriormente, la actividad regulatoria del SNP principal rs6061231 estuvo dominada por rs67941642, y el alelo rs6061231-C confirió una actividad significativamente mayor que el alelo A.

Los SNPs identificados se asocian con la metástasis del CRC

Para investigar más a fondo la relevancia clínica de los SNPs rs606

Utilizando el sistema SNP-STARR-seq, realizamos un cribado de SNPs que presentan diferencias alelo-específicas en la actividad de los potenciadores entre las líneas celulares SW480 y SW620 (tabla S3). Los paSNPs identificados en las células SW620 fueron predominantemente variantes funcionales, con un 88,9% respaldado por evidencia de eQTL y un 72,3% que se predijo que residían dentro de sitios de unión a factores de transcripción (TF) (fig. S3A). La comparación cruzada de los resultados de las dos líneas celulares identificó 3136 SNPs que exhibieron actividad de potenciador alelo-específica específicamente en las células SW620 (fig. 4A y tabla S6). El análisis HOMER predijo un enriquecimiento significativo de motivos de unión para los TFs de la familia bZIP alrededor de estos SNPs (fig. 4, B y C), una familia estructural de TFs predominantemente asociada con la metástasis del CRC (^[16], ^[22]–^[25]).

Luego, buscamos refinar aún más el mapeo de variantes genéticas funcionales y loci que contribuyen a la metástasis del CRC. El análisis de RNA-seq que comparó las células SW480 y SW620 identificó 2560 genes que se sobreexpresaron de manera significativa y específica en SW620, lo que sugiere su fuerte asociación con la metástasis del CRC [P ajustada < 0,05, cambio de pliegue logarítmico (log2FC) > 0,75; fig. S3B y tabla S7]. Los análisis de enriquecimiento funcional [Ontología Genética (GO) y Enciclopedia de Genes y Genomas de Kioto (KEGG)] confirmaron que estos genes sobreexpresados estaban estrechamente relacionados con procesos relacionados con la metástasis, como la migración celular y las vías de señalización de Wnt (fig. S3C). Entre los 3136 SNPs reguladores específicos de SW620, 2517 (80,3%) se ubicaron dentro de ±1 Mbp de los sitios de inicio de la transcripción de estos genes sobreexpresados en relación con la metástasis, lo que sugiere una relación regulatoria entre las variantes genéticas y estos genes. Al integrar estos hallazgos con los picos de ChIP-seq de H3K27ac específicos de SW620, demostramos que 1697 de estos SNPs se situaron en potenciadores activados por la metástasis. En conjunto, estos resultados subrayan que SNP-STARR-seq es una herramienta eficaz para identificar variantes genéticas funcionales de manera específica para cada tipo de célula, lo que facilita el descubrimiento de SNPs de riesgo asociados con la metástasis del CRC.

El papel del SNP rs1962004 en la regulación de la expresión de LRRC61

Dado el marcado enriquecimiento de motivos de unión a TFs de la familia bZIP en estos loci, hipotetizamos que los SNPs que interrumpen la unión a bZIP impulsan el secuestro de potenciadores durante la progresión metastásica. Por lo tanto, priorizamos 22 SNPs asociados con la metástasis de alta confianza (fig. 4D y tabla S6) basándonos en cuatro criterios, que incluyen: (i) actividad de potenciador alelo-específica significativa en SW620 (pero negativa en SW480), (ii) ubicación dentro de los picos de H3K27ac específicos de SW620, (iii) proximidad (±1 Mbp) a los genes sobreexpresados en relación con la metástasis y (iv) capacidad de interrumpir los supuestos sitios de unión a TFs de la familia bZIP. Entre estas variantes priorizadas, encontramos que un SNP, rs1962004, presentaba diferentes genotipos entre las células SW480 y SW620 (SW480, G/G; y SW620, G/A; fig. 4E), lo que indica que el alelo "A" específico de SW620 era muy probablemente una mutación adquirida somáticamente y podría contribuir a la metástasis de las células de CRC. Por lo tanto, seleccionamos el SNP para una posterior caracterización mecanicista. En comparación con las células SW480, las células SW620 mostraron un enriquecimiento adquirido de H3K27ac en este locus, lo que revela la actividad del potenciador de la vecindad. Mientras tanto, observamos una expresión significativamente más fuerte del gen adyacente LRRC61 en las células SW620 en comparación con las células SW480, lo que indica la posible relación regulatoria entre el potenciador específico de la metástasis y la expresión del gen LRRC61 (fig. 4F).

Como se discutió anteriormente, las mutaciones no codificantes generalmente conducen a una unión alterada de los TFs y a una expresión aberrante del gen diana. Se predijo que el SNP rs1962004 residía dentro de un sitio de unión putativo de los TFs de la familia bZIP, JUN y JUND. Por lo tanto, primero confirmamos la unión significativa de JUN y JUND a este locus utilizando ChIP-qPCR (fig. 4G). La redundancia funcional se alinea con los mecanismos establecidos de los TFs bZIP, donde JUN y JUND forman con frecuencia heterodímeros como parte del complejo AP-1 para regular cooperativamente los genes diana involucrados en la oncogénesis y la metástasis (^[26]–^[28]). Cabe destacar que estos TFs prefirieron el alelo "A" adquirido por SW620, lo que sugiere que se creó un nuevo sitio de unión de JUN y JUND en el locus durante la metástasis cuando SW620 experimentó una transición de "G" a "A" de rs1962004. Esto también fue consistente con los cambios epigenéticos observados en el locus circundante de las células metastásicas de SW480 a SW620 (fig. 4F). Además, un ensayo de reportero de luciferasa demostró una actividad significativamente mayor del fragmento genómico que contiene el alelo "A" en comparación con el alelo "G" (P < 0,001, fig. 4H).

Además, utilizando la edición del genoma basada en CRISPR-Cas9, generamos dos clones celulares con genotipos homozigotos del SNP rs1962004: G/G y A/A, respectivamente. Nuestros resultados mostraron de manera consistente que las células con el genotipo A/A exhibieron una capacidad proliferativa y migratoria significativamente mayor en comparación con las células SW620 con el genotipo heterocigoto rs1962004 (es decir, G/A), que, a su vez, exhibieron una mayor capacidad pro-metastásica que las células homozigotas G/G (fig. 4, I a K). Complementariamente, la deleción de la región del potenciador de 104 pares de bases que abarca rs1962004 (fig. S3, D y E) atenuó significativamente la proliferación en las células SW620. Mientras tanto, la expresión de LRRC61 se redujo significativamente en las células con el potenciador eliminado en comparación con los controles de tipo salvaje (P < 0,01, fig. S3F). En conjunto, estos hallazgos demuestran que tanto el alelo "A" adquirido como el elemento del potenciador en sí impulsan mecánicamente la agresividad metastásica. La unión de JUN o JUND, causada por la transición de "G" a "A" en rs1962004, mejoró posteriormente la actividad cis-reguladora y aumentó la expresión de LRRC61.

El papel metastásico de LRRC61

Para establecer aún más la relación regulatoria entre el SNP y LRRC61, consultamos los datos de eQTL del conjunto de datos TCGA-CRC y revelamos que LRRC61 exhibió una expresión correlacionada con los genotipos de rs1962004 (P = 7,19 × 10^19). Es decir, los individuos que portaban el genotipo A/A tenían una expresión significativamente mayor de LRRC61, seguido por los portadores de A/G y G/G (fig. 5A). En consonancia con las tendencias fenotípicas más amplias observadas a nivel de la población, las células SW620 editadas con un genotipo A/A también mostraron niveles elevados de LRRC61 en comparación con las células de la misma procedencia que portaban el genotipo G/G (fig. 5B). El análisis previo reveló un enriquecimiento relacionado con la metástasis de JUN y JUND en el locus rs1962004. Aquí, realizamos la supresión mediada por siRNA de JUN y JUND, ambos de los cuales atenuaron la expresión de LRRC61 (fig. 5, C y D).

La superfamilia de proteínas que contienen repeticiones ricas en leucina (LRRC) tiene funciones complejas y diversas en varios tumores sólidos (como el sarcoma, el melanoma y el glioblastoma), y algunos de sus miembros (como LRRC15 y LRRC32) impulsan la invasión y la metástasis del tumor al regular la señalización del factor de crecimiento transformante-β, promover la activación de fibroblastos asociados con el cáncer e influir en el microambiente inmunitario (^[29], ^[30]). Sin embargo, el papel de LRRC61 en el cáncer, especialmente en el CRC, no se había comprendido bien. Aquí, pretendimos investigar su posible participación en el impulso de los fenotipos metastásicos. La supresión de LRRC61 (P < 0,001) utilizando ARN de interferencia corta (shRNA) condujo a una disminución de la capacidad de formación de colonias y la migración (fig. 5, E a G), que son fenotipos estrechamente asociados con la metástasis. La metástasis causa significativamente la mortalidad del paciente y un mal pronóstico. El seguimiento clínico de los pacientes con CRC reveló que la expresión de LRRC61 era mayor en los tejidos tumorales en comparación con los tejidos normales (fig. 5H), y los pacientes con una expresión elevada de LRRC61 generalmente tenían un peor pronóstico (fig. 5I). Un modelo genético para esta regulación propone que el alelo de riesgo rs1962004-A, al mejorar la unión de JUN/JUND y aumentar la expresión de LRRC61, subyace al aumento del potencial metastásico y al peor pronóstico del paciente (fig. 5J).

Methyl-STARR-seq identifica elementos potenciadores dependientes de la metilación

Además de las variantes genéticas que influyen en los fenotipos de la enfermedad, las variaciones epigenéticas también se reconocen ampliamente como factores etiológicos clave que afectan la salud humana. La metilación del ADN representa un mecanismo epigenético crítico que regula la expresión génica. Si bien la investigación existente se centra predominantemente en la metilación del promotor, los roles regulatorios y el mapeo genómico preciso de la metilación dentro de las regiones de los potenciadores, que constituyen una proporción sustancial del epigenoma, siguen estando insuficientemente caracterizados. Para abordar esta brecha, diseñamos una biblioteca de oligonucleótidos dirigida a regiones de potenciadores metiladas diferencialmente en el CRC y ventanas de 120 nt centradas en los picos de H3K27ac. Esta biblioteca se sintetizó y clonó en el vector CpG-free Methyl-STARR-seq (^[14]). Tras el tratamiento enzimático in vitro con glicerol como control o M.SssI, una enzima que puede catalizar la metilación completa de la secuencia CpG únicamente en los fragmentos insertados (^[31]), generamos versiones no metiladas y completamente metiladas de los mismos grupos de plásmidos. La tasa de conversión y el estado de metilación se confirmaron mediante la secuenciación de bisulfito (BS-seq; fig. S4, A y B). Estos constructos se transfectaron en las líneas celulares HCT116, SW480 y SW620, respectivamente, para cuantificar su actividad de potenciador. Posteriormente, se construyeron bibliotecas de ARNm y ADN de entrada y se sometieron a secuenciación y análisis bioinformáticos (esquema del flujo de trabajo experimental en la fig. 1). Al calcular la actividad diferencial del potenciador (EAS) entre los estados metilados y no metilados, identificamos los elementos reguladores más significativamente influenciados por la metilación del ADN en cada línea celular (HCT116, n = 1763; SW480, n = 3588; y SW620, n = 4112; fig. S4, C y D, y tabla S4).

Tras la metilación, aproximadamente el 60% de los fragmentos probados exhibieron una actividad de potenciador reducida en las células HCT116, en comparación con aproximadamente el 40% en las células SW480 y SW620, lo que destaca la heterogeneidad específica de la línea celular (fig. 6A). Para los elementos comunes a HCT116 y SW480, la validación utilizando un vector de reportero de Lucia CpG-free (InvivoGen) demostró una alta concordancia con los resultados de la secuenciación de alto rendimiento (fig. 6B y fig. S4E). Las alteraciones dependientes de la metilación en la actividad del potenciador observadas probablemente se debieron a la sensibilidad a la metilación de los TFs. El análisis del potencial de unión a TFs reveló sitios de unión predichos dentro de más del 95% de estos elementos (tabla S8). Cabe destacar que casi el 80% de estos sitios contenían motivos ricos en CpG dentro de sus secuencias de unión a TFs predichas, lo que sugiere que la metilación de CpG podría afectar directamente la ocupación de los TFs (fig. 6C). Además, los TFs enriquecidos en estos elementos reguladores identificados por Methyl-STARR-seq mostraron una fuerte concordancia con las categorías de TFs methylplus (definidas como la metilación de CpG dentro de los sitios de unión aumenta la afinidad de unión) y methylminus (definidas como la metilación dentro de los sitios de unión de CpG disminuye la afinidad de unión) definidas por Methylation Systematic Evolution of Ligands by Exponential Enrichment (Methyl-SELEX) (^[32]), lo que confirma la fiabilidad de nuestro sistema de cribado y también el papel de los TFs en la modulación de la actividad del potenciador dependiente de la metilación (fig. 6D). Realizamos además un análisis de enriquecimiento de motivos HOMER completo para los elementos methylminus frente a methylplus identificados en las líneas celulares HCT116, SW480 y SW620, así como específicamente para los 3008 elementos reguladores sensibles a la metilación adquiridos por la metástasis únicos de SW620. Nuestros resultados revelaron que los elementos methylminus en las líneas celulares HCT116 y SW620 más agresivas se enriquecieron preferentemente con TFs de la familia ETS, que son impulsores bien establecidos de la invasión y la metástasis (^[33], ^[34]). En contraste, la línea primaria SW480 mostró un mayor enriquecimiento de los motivos bZIP/AP-1. Este notable patrón específico del contexto, en particular el fuerte enriquecimiento del factor ETS en los elementos sensibles a la metilación adquiridos por la metástasis de SW620, proporciona un vínculo mecanicista convincente entre la reprogramación epigenética durante la progresión metastásica y la activación de programas transcripcionales pro-metastásicos (fig. S5A).

Posteriormente, nos centramos en las DMR que podrían desempeñar un papel fundamental en la patogénesis del CCRC. Integrando el análisis de motivos TF con los datos de modificación de H3K4me1, identificamos seis sitios CpG funcionalmente putativos, incluidos cg25982657, cg08640619, cg23516310, cg14302214, cg14179401 y cg14106234, ninguno de los cuales había sido caracterizado previamente (fig. 6E y tabla S9). El análisis de los datos de metilación de muestras tumorales y normales emparejadas de 14 tipos de cáncer reveló que cg25982657 y cg08640619 estaban hipermetilados específicamente en el CCRC (| | > 0,25, tasa de descubrimiento falso (FDR) < 0,05; fig. S5B y tabla S10].

Utilizando cg08640619 como un sitio representativo, realizamos una validación funcional. Este sitio CpG se encontraba en el primer intrón del gen KIRREL1, que exhibía una actividad potenciadora significativa (fig. 6F). La predicción in silico indicó una posible unión de RUNX2, un TF que se sabe que se une preferentemente al ADN no metilado. Para confirmar el papel de la metilación del ADN en la unión de RUNX2 a este locus, utilizamos la metilación dirigida mediada por CRISPR-Cas9 inactivo (dCas9) en las células HCT116 y SW480. Posteriormente, el ChIP-qPCR demostró una reducción significativa del enriquecimiento de RUNX2 en el locus hipermetilado en las células HCT116 (P < 0,001, fig. 6G). La inhibición de la expresión de RUNX2 provocó una regulación a la baja significativa de sus genes adyacentes, KIRREL1 y ETV3, dentro del mismo TAD, lo que confirmó que RUNX2 se unía preferentemente a la secuencia no metilada y regulaba positivamente KIRREL1 y ETV3 (fig. 6H). Tras observar una correlación negativa significativa entre la metilación de cg08640619 y la expresión de KIRREL1/ETV3 en muestras clínicas (fig. 6I), buscamos validar funcionalmente este hallazgo. Utilizando la hipermetilación dirigida mediada por CRISPR-dCas9-DNMT3A, confirmamos que el aumento de la metilación provoca directamente una regulación a la baja significativa de ambos genes, lo que demuestra una relación regulatoria causal (fig. S5C). Se descubrió que la expresión de KIRREL1 estaba significativamente aumentada en los tejidos tumorales en comparación con los tejidos normales correspondientes en las muestras emparejadas de adenocarcinoma de colon (COAD) del TCGA (fig. S5D). Este aumento de la expresión puede estar asociado con el estado de hipometilación de cg08640619 en los tumores de CCRC. El análisis pronóstico de los datos clínicos reveló que la hipometilación de cg08640619 está significativamente asociada con un mal pronóstico del paciente en la enfermedad en estadio avanzado (estadios III a IV) (fig. 6J).

Dado los patrones de hipometilación específicos del cáncer (fig. S5B), hipotetizamos que cg08640619 y cg25982657 podrían servir como biomarcadores sensibles para la detección temprana del CCRC. Para evaluar el potencial diagnóstico de los sitios CpG identificados, aplicamos tanto la regresión logística (LogR) como los clasificadores de máquinas de vectores de soporte (SVM) a muestras de 14 tipos de cáncer diferentes. Además, se realizó una validación cruzada de 10 pliegues para garantizar la solidez del modelo (fig. S5, E, G y H). Los resultados mostraron que ambos sitios CpG discriminaban eficazmente los tumores de CCRC de los tejidos normales, pero no tan bien para otros tumores (fig. 6K y fig. S5F). Los variantes epigenéticos funcionales detectados por Methyl-STARR-seq estaban altamente correlacionados con el riesgo de enfermedad, lo que les confiere un gran potencial clínico para el diagnóstico específico y sensible del cáncer.

Methyl-STARR-seq identifica la metilación dependiente de la metástasis

De manera similar a nuestro cribado de SNPs funcionales asociados con la metástasis, identificamos 3008 regiones reguladoras sensibles a la metilación que se desregulan específicamente en las células SW620 (fig. 6L y tabla S11). El análisis integrativo de los perfiles de modificación de H3K27ac y los datos de RNA-seq en ambas líneas celulares identificó 182 fragmentos genómicos que exhibían una actividad reguladora mejorada en el estado de hipometilación (fig. 6M). Estos loci mostraron simultáneamente un enriquecimiento significativamente elevado de H3K27ac en las células SW620 en comparación con las células SW480 y se ubicaron dentro de ± 1 Mb de genes que se expresaban significativamente más en esta línea metastásica.

Tomemos un sitio putativo, cg02205193, como ilustración. Estos dinucleótidos CpG se encuentran cerca del promotor de LINC00460, un gen que codifica un ARN endógeno competidor que se ha informado que promueve la metástasis del CCRC en varios estudios (^[35]–^[37]). La región genómica que encierra cg02205193 exhibió señales de H3K27ac significativamente más altas en las células SW620 en comparación con las células SW480 (fig. 6N). En consecuencia, LINC00460 se expresa significativamente más en los pacientes con CCRC metastásico, como se observó en la línea celular metastásica SW620, lo que podría atribuirse a su estado de hipometilación (^[38]). Además, la metilación en cg02205193 influyó supuestamente en la unión de GFI1, un TF que no se une a la metilación, según los datos de Methyl-SELEX (^[32]). Específicamente, el estado de hipometilación de cg02205193 facilita una mayor captación de GFI1, lo que, a su vez, activa la expresión del gen diana LINC00460 y promueve la metástasis del CCRC. Sobre la base de Hi-C y la expresión diferencial de genes entre las células SW480/SW620, también identificamos regiones reguladoras que encierran sitios CpG previamente no caracterizados cuyo estado de metilación se asoció con la expresión de los genes TCF12 y HNF1, dos genes establecidos implicados en la metástasis del CCRC (fig. S5I) (^[39], ^[40]).

Por lo tanto, al integrar el cribado funcional de Methyl-STARR-seq con la validación multiómica, pudimos identificar eficazmente posibles biomarcadores diagnósticos. Este enfoque no solo delimita un paisaje regulatorio preciso de la desregulación epigenética en el CCRC, sino que también proporciona una base mecánica para la identificación de biomarcadores diagnósticos y terapéuticos previamente no identificados con una amplia aplicabilidad clínica.

DISCUSIÓN

En este estudio, establecimos una vía integrativa de genómica funcional para priorizar sistemáticamente las variantes genéticas y epigenéticas no codificantes que impulsan la patogénesis y la metástasis del CCRC. Al combinar SNP-STARR-seq (cribado de 30.790 SNPs) y Methyl-STARR-seq (perfilado de >134.000 sitios CpG) en células primarias (HCT116 y SW480) y metastásicas (SW620) de CCRC, mapeamos las actividades potenciadoras dependientes de alelos y de metilación a una escala sin precedentes. Identificamos 1969 SNPs funcionales en HCT116 y 1846 en las células primarias de CCRC SW480, con 922 variantes comunes que impulsan la desregulación del potenciador. En particular, se descubrieron 3136 SNPs adquiridos en la metástasis y 3008 elementos sensibles a la metilación en las células metastásicas SW620. La validación mecánica confirmó que TAF4 es un supresor tumoral regulado por rs67941642/rs6061231. Encontramos que una variante adquirida en la metástasis, rs1962004 (alelo A), mejora la expresión de LRRC61 a través de la unión de JUN/JUND, lo que promueve la progresión del CCRC y se correlaciona con un mal pronóstico. Además, identificamos loci hipermetilados, como cg08640619 y cg25982657, como biomarcadores diagnósticos de alto rendimiento (AUROC > 0,96). Este enfoque integrativo une tanto la variación genética como la epigenética a la reprogramación transcripcional en la patogénesis y la metástasis del CCRC.

Nuestro estudio, junto con trabajos recientes en células normales primarias (^[10]), destaca que la lectura funcional de las variantes no codificantes está profundamente influenciada por el contexto celular. La actividad reguladora mejorada que observamos en las líneas celulares de CCRC probablemente refleja el entorno transcripcional y epigenético aberrante de los tumores avanzados, lo que revela un espectro más amplio de variantes funcionalmente relevantes. Esta especificidad del contexto, junto con nuestro modelo metastásico isogénico, fue fundamental para identificar los impulsores reguladores de la metástasis. La caracterización funcional de los SNPs no codificantes en el cáncer ha sido revolucionada por las tecnologías clave, MPRA o STARR-seq. MPRA, pionera como un "detector de mentiras reguladoras", permite la prueba simultánea de miles de fragmentos de ADN para la actividad del potenciador al acoplar secuencias candidatas a los reporteros con código de barras (^[41], ^[42]). Si bien estos métodos han vinculado con éxito los loci de riesgo de GWAS a los mecanismos reguladores en varios contextos (^[42]–^[46]), su aplicación en el CCRC, particularmente con respecto a la metástasis, ha sido limitada y restringida por las opciones de diseño específicas. Por ejemplo, el estudio de SNPs de CCRC más completo hasta la fecha utilizó MPRA, pero se centró exclusivamente en las variantes dentro de los loci de riesgo de GWAS mapeados estadísticamente, probando solo 8880 variantes en las células primarias de CCRC (^[46]). Este enfoque, si bien es valioso, inherentemente pasó por alto el vasto panorama de los SNPs reguladores funcionales que se encuentran fuera de las regiones implicadas en GWAS y el contexto crítico de la progresión metastásica. Nuestro estudio avanza fundamentalmente más allá de estas limitaciones a través de una estrategia conceptualmente distinta y más amplia: (i) Cribado de variantes a doble escala: implementamos una estrategia de cribado de SNPs-STARR-seq a doble escala. A diferencia de los estudios anteriores que priorizaban solo las variantes vinculadas a GWAS, interrogamos sistemáticamente tanto los SNPs asociados a GWAS (6363) como los SNPs localizados en potenciadores independientes de los éxitos de GWAS (24.526), para un total de 30.790 variantes únicas. Este enfoque imparcial amplió significativamente el alcance más allá de los loci de riesgo conocidos, capturando una imagen más completa de la variación reguladora no codificante funcional dentro de los potenciadores del CCRC. (ii) Incorporación del contexto metastásico: en particular, realizamos el cribado no solo en las líneas celulares primarias de CCRC (HCT116 y SW480), sino también en la línea metastásica emparejada isogénicamente SW620, derivada de la metástasis del mismo paciente en el ganglio linfático. Este modelo emparejado único es fundamental, ya que nos permitió comparar directamente los paisajes de actividad del potenciador y señalar 3136 SNPs reguladores adquiridos en la metástasis que son funcionalmente específicos del nicho metastásico, una dimensión completamente ausente en los estudios de genómica funcional a gran escala anteriores del CCRC. La caracterización funcional de estas variantes específicas de la metástasis reveló un enriquecimiento notable y significativo de los motivos de unión de los TFs de la familia bZIP (por ejemplo, JUN y JUND). Este hallazgo implica fuertemente que la desregulación de los complejos bZIP/AP-1, los principales reguladores de la transición epitelial-mesenquimal y la metástasis (^[16], ^[22]–^[25]), es una consecuencia clave de la evolución somática y reguladora durante la propagación metastásica del CCRC. La identificación de este enriquecimiento específico de la familia de TFs entre las variantes adquiridas en la metástasis es una idea original que se habilita de forma única mediante nuestro cribado en el contexto metastásico.

Se ha establecido que la metilación del ADN es un mecanismo epigenético crucial en el diagnóstico del cáncer (^[47]–^[49]). Sin embargo, la investigación existente sobre la metilación del cáncer se ha centrado principalmente en las regiones promotoras, pasando por alto el potencial significativo de las regiones potenciadoras, que representan una proporción sustancial del epigenoma (^[9]). En este estudio, utilizamos Methyl-STARR-seq (^[14]), que se ha aplicado con menos frecuencia en el contexto del cáncer, para investigar sistemáticamente la metilación de la región potenciadora en el CCRC. En particular, dado que tanto los SNPs como los sitios de metilación se seleccionaron a partir de las mismas regiones potenciadoras, este enfoque integrado proporciona información mecánica sobre la arquitectura reguladora de los potenciadores del CCRC, lo que delimita los segmentos genómicos influenciados por la variación genética frente a aquellos que se modulan preferentemente mediante alteraciones epigenéticas (tabla S1). Al integrar el enfoque de Methyl-STARR-seq con los clasificadores de aprendizaje automático, identificamos dos posibles biomarcadores de alerta temprana, cg08640619 y cg25982657. Estas dos sondas metiladas diferencialmente representan las variaciones de metilación específicas del CCRC que exhiben un notable potencial diagnóstico (fig. 6K, fig. S5F y tabla S10).

Nuestro marco integrado de genómica funcional establece un modelo para descifrar los factores no codificantes en diversos tipos de cáncer. La metodología y las conclusiones derivadas de este estudio centrado en el CCRC ofrecen implicaciones clave para la traslación a otros tipos de cáncer: aplicabilidad pan-cáncer de la estrategia de cribado. La plataforma dual SNP-/Metil-STARR-seq, junto con la caracterización multi-ómica, es fácilmente adaptable a otros cánceres epiteliales. Por ejemplo, la selección de variantes centrada en potenciadores (priorizando las regiones definidas por H3K27ac/ATAC) evita las limitaciones del GWAS en cánceres poco estudiados (por ejemplo, adenocarcinoma pancreático), lo que permite el descubrimiento de novo de factores reguladores. La manipulación in vitro de la metilación mediante M.SssI permite la elaboración sistemática de perfiles de potenciadores sensibles a la metilación en cánceres con una importante desregulación epigenética [por ejemplo, subtipos de glioblastoma con mutaciones en la isocitrato deshidrogenasa (IDH)].

Desarrollo de biomarcadores predictivos y diagnósticos: Nuestra estrategia de descubrimiento de biomarcadores, desde el cribado funcional hasta la detección de SNPs causales o regiones de metilación diferencial de potenciadores (eDMR), proporciona una plantilla para otros tipos de cáncer: los potenciadores hipermetilados específicos de tejido (por ejemplo, cg08640619 en nuestro estudio) podrían utilizarse en cánceres que carecen de biomarcadores sensibles (por ejemplo, cáncer de páncreas en estadio temprano). El rendimiento diagnóstico excepcional (AUROC > 0,96) logrado por nuestros marcadores de metilación de potenciadores específicos de CCRC, cg08640619 y cg25982657, supera significativamente el rango de rendimiento típico que se informa para los biomarcadores de metilación basados en promotores en la detección temprana de CCRC, que a menudo presentan valores de AUROC entre 0,80 y 0,90 (^[50]–^[52]). Esta mejora sustancial subraya el poder del cribado funcional de variantes epigenéticas dentro de los elementos potenciadores, en lugar de depender únicamente de la elaboración de perfiles centrados en el promotor. Estas firmas epigenómicas informadas funcionalmente no solo proporcionan una mayor precisión diagnóstica, sino que también identifican elementos reguladores vinculados mecánicamente a la patogénesis de la enfermedad.

Si bien nuestro enfoque integrado proporciona una anotación funcional completa de las variantes no codificantes en el CCRC, varias limitaciones deben tenerse en cuenta. En cuanto a las limitaciones del diseño de la biblioteca impuestas por las limitaciones de la síntesis de oligonucleótidos, los fragmentos genómicos capturados (excluidas las secuencias de los adaptadores) se limitaron a ventanas de 120 pb. Esta escala de 120 pb puede ser insuficiente para recapitular completamente las interacciones de la cromatina a larga distancia o los efectos del dominio topológico. Además, para la elaboración de perfiles de metilación, si bien las 18.880 regiones de potenciadores albergaban más de 950.000 sitios CpG, la variación sustancial en las longitudes de los potenciadores (que oscilan entre cientos y miles de pares de bases) requirió que nos centráramos en segmentos de 120 pb sintetizados centrados en los picos de H3K27ac. Reconocemos que nuestro diseño dirigido necesariamente omite ciertas clases de loci funcionales: potenciadores latentes o en reposo: las regiones que no están marcadas por H3K27ac en las líneas celulares que analizamos (por ejemplo, potenciadores débiles y potenciadores activos en otros subtipos de CCRC o bajo estímulos específicos) estarían ausentes de nuestra biblioteca. Además, nuestros modelos celulares, si bien proporcionan condiciones experimentales controladas, carecen de la complejidad de los microambientes tumorales (por ejemplo, las interacciones inmuno-estromales) que pueden modular la actividad de los potenciadores in vivo. Los estudios futuros deben priorizar: (i) estrategias de fragmentación in situ para construir bibliotecas de tamaño diverso que abarquen la arquitectura completa de los potenciadores; (ii) enfoques multi-ómicos de célula única para analizar la heterogeneidad regulatoria específica de alelos y metilación durante la progresión metastásica; y (iii) si bien hemos validado el potencial diagnóstico de cg08640619 y cg25982657 en muestras de tejido, su aplicación como biomarcador en biopsias líquidas requiere una mayor investigación y desarrollo. Abordar estas dimensiones permitirá salvar la brecha entre los mapas de potenciadores y los contextos terapéuticos aplicables.

MATERIALES Y MÉTODOS

Diseño de oligos

Establecimos un conjunto integral de SNPs a través de criterios integrados de funcionalidad y relevancia para la enfermedad: (i) aplicamos un filtrado estricto de la frecuencia del alelo menor (MAF) [0,05 en AMR/EUR/EAS/SAS en el Proyecto de los 1000 Genomas (^[15])]; (ii) ampliamos 238 variantes de GWAS de CCRC mediante la imputación de proxy de desequilibrio de enlace (LD) (R2 ≥ 0,5, ±500 kb), lo que dio como resultado 6125 SNPs adicionales; y (iii) priorizamos de forma independiente 24.526 SNPs localizados en potenciadores definidos por la coaccesibilidad de la cromatina [intersección de H3K27ac ChIP-seq (ENCSR661KMA/ENCSR000EUT) y picos de ATAC-seq (ENCSR872WGW) en células HCT-116]. El conjunto de unión comprendió 30.790 SNPs únicos (6363 derivados de GWAS y 24.526 basados en potenciadores, con 99 superposiciones).

Para los oligos específicos de metilación, seleccionamos (i) los CpGs del núcleo del potenciador: más de 134.000 dinucleótidos CpG dentro de ±60 pb de los picos de H3K27ac; y (ii) CpGs diferenciales de CCRC (| diferencia | > 0,3, FDR < 0,05): 908 sondas localizadas en potenciadores de las matrices Illumina 450K del conjunto TCGA-CRC, GSE139404, GSE77955, GSE53051, GSE48684, GSE42752 y GSE40055. El diseño de los oligonucleótidos se adaptó a un sistema de secuenciación de Beijing Genomics Institute (BGI), y las dos bibliotecas se sintetizaron por separado por CustomArray. En resumen, cada oligonucleótido contiene 120 pb de secuencia genómica que encierra los sitios de SNP o CpG y secuencias de flanqueo constantes (en el extremo 5', 5'-ATGGCTACGATCCGACTT; y en el extremo 3', AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA-3'), que se utilizaron para el clonado.

Preparación del plásmido SNP-STARR-seq

El vector hORI-STARR-seq fue proporcionado por A. Stark (Instituto de Investigación de Patología Molecular, Viena, Austria). La biblioteca de SNPs amplificada por PCR (20 ng) se insertó en los sitios Age I/Sal I del plásmido (100 ng) utilizando un clonado basado en recombinación altamente eficiente (ClonExpress, Vazyme, China). Para evitar sesgos durante el clonado, realizamos un total de 20 reacciones de recombinación separadas y agrupamos cada cuatro reacciones. Después de una purificación de ADN con cuentas AMPure XP, utilizamos los vectores que contenían los SNPs (2 μl) para transformar Escherichia coli DH5α extra-competente (DH5α; 50 μl). Nuevamente, realizamos 20 reacciones de transformación separadas y agrupamos cada cuatro transformaciones. Luego, los plásmidos se extrajeron y purificaron a una alta concentración (>1 μg/μl). Se realizó una amplificación por PCR a partir de los plásmidos y se secuenciaron para 2 × 100 ciclos de lectura emparejada como control de entrada.

Preparación del plásmido Metil-STARR-seq

El pmSTARRseq1 fue un regalo de J. Tung. Los métodos de clonado son similares a los del SNP-STARR-seq mencionado anteriormente, excepto que la biblioteca de metilación amplificada por PCR (10 ng) se insertó en los sitios Spe I/Nco I del plásmido (100 ng). Los plásmidos se transformaron en células GT115 extra-competentes (con selección con zeocina a 50 μg/ml). Luego, metilamos la mitad de la biblioteca de plásmidos con M.SssI [ADN, 4 U/μg; New England Biolabs (NEB)] y tratamos la otra mitad como control (utilizando el mismo protocolo de reacción, pero reemplazando la enzima con agua destilada). Para confirmar el estado de metilación del ADN del plásmido, primero digerimos el ADN del plásmido (1 μg) incubándolo (a 37 °C) con Xba I y Nhe I en un total de cinco reacciones. Aproximadamente 1 ng de fragmento de proteína fluorescente verde no metilada se añadió a cada muestra para evaluar la eficiencia de la conversión con bisulfito. Realizamos reparación de extremos, adición de cola A y ligación de adaptadores. Luego, las muestras se sometieron a conversión de sodio bisulfito utilizando un kit EZ DNA Methylation-Gold (Zymo Research). Las bibliotecas se amplificaron por PCR (ADN polimerasa KAPA HiFi Uracil+) y se secuenciaron en la plataforma Illumina para 2 × 100 ciclos de lectura emparejada.

Cultivo celular y transfección

Cada grupo de oligos se transfectó en células embrionarias humanas (HEK 293T) o líneas celulares de CCRC (HCT116, SW480 y SW620) utilizando jetOPTIMUS (Polyplus). Las células se cultivaron en condiciones normales con 5 % de CO2 a 37 °C. Específicamente, se mezclaron 2 μg de plásmidos con 3 μl de reactivos de transfección para la transfección en células cultivadas en un pozo de una placa de seis pocillos. Agrupamos cada seis pocillos (~2 × 107 células) para una repetición biológica.

Extracción de ARNm y secuenciación

Treinta y seis horas después de la transfección, se extrajo el ARN total con un kit RNeasy (TAKARA) y se enriqueció el ARNm con un módulo de aislamiento magnético de ARNm poli(A) (NEB, E7490L). El ARNm residual se eliminó mediante tratamiento con ribonucleasa (RNasa) A/H. El ADNc se purificó y recuperó con cuentas AMPure XP y se ligaron adaptadores. Luego, el ADN se amplificó con cebadores de secuenciación BGI (BGI-system F y BGI-Barcode R), utilizando diferentes cebadores posteriores para diferentes muestras para permitir el multiplexado de muestras. Después de la amplificación, el ADN se enriqueció aún más con cuentas AMPure XP antes de agruparlo para la secuenciación.

Ensayos de reportero de luciferasa dual

Para investigar el impacto funcional de los SNPs o los sitios de metilación en la actividad de los potenciadores, se realizaron ensayos de reportero de luciferasa dual. Brevemente, fragmentos de ADN de 120 pb que contenían diferentes alelos de SNPs o fragmentos que contenían CpG se clonaron en el vector de reportero de luciferasa pGL4.23 (sitios Kpn I/Xho I) o en el vector pCpG-free-lucia (sitios BamH I/Spe I), respectivamente. El plásmido recombinante pGL4.23 o pCpG-free-lucia (reportero de luciferasa de luciérnaga) se cotransfectó con el vector pGL4.74 para normalizar las variaciones en la eficiencia de la transfección y la viabilidad celular. Después de 48 horas de transfección, se midió la actividad de la luciferasa utilizando el kit de ensayo de luciferasa Dual-Lumi (Beyotime, RG088S). Actividad relativa del potenciador = Luminiscencia de luciérnaga/Luminiscencia de renilla. Los cebadores de clonado, así como las secuencias de los sitios variantes detectados, se enumeran en la tabla S13. Los datos originales de la luminiscencia de luciérnaga y renilla se enumeran en la tabla S14.

Secuenciación de ARN

Se extrajo el ARN total de las células utilizando el kit de extracción de ARN universal MiniBEST de TAKARA. Se enriqueció el ARN poliadenilado a partir de 5 μg de ARN total, seguido de la elución directa utilizando un sistema de transcripción inversa preformulado para la síntesis de ADNc. El ARN residual se eliminó mediante tratamiento con ribonucleasa (RNasa) A/H. El ADNc se purificó y recuperó con cuentas AMPure XP y se ligaron adaptadores. Luego, el ADN se amplificó con cebadores de secuenciación BGI (BGI-system F y BGI-Barcode R), utilizando diferentes cebadores posteriores para diferentes muestras para permitir el multiplexado de muestras. Después de la amplificación, el ADN se enriqueció aún más con cuentas AMPure XP antes de agruparlo para la secuenciación.

Genotipado

Se extrajo el ADN genómico de las células siguiendo el protocolo del fabricante del kit de extracción de ADN. Se diseñaron pares de cebadores que flanquean el locus de SNP objetivo (tabla S15) y se utilizaron para la amplificación por PCR con ADN polimerasa Taq. Los productos de PCR se enviaron directamente a una empresa de biotecnología para la secuenciación de Sanger. El genotipo del locus de SNP se determinó analizando los picos del cromatograma de secuenciación en la posición objetivo (SnapGene).

Edición genética mediada por CRISPR-Cas9

Eliminación

Para dilucidar el papel funcional de regiones específicas que contienen SNPs como posibles potenciadores, utilizamos la edición genética mediada por CRISPR-Cas9 para generar la eliminación dirigida del fragmento del sitio de SNP. Este enfoque permitió evaluar la actividad del potenciador mediante el análisis cuantitativo de los cambios en la expresión génica aguas abajo.

Se diseñaron dos ARN guía únicos (sgARN) que flanquean la región del polimorfismo de un solo nucleótido (SNP) objetivo, utilizando la plataforma de diseño CRISPR CHOPCHOP (https://chopchop.cbu.uib.no/). Los sgARN se seleccionaron en base a criterios estrictos, que incluyen una alta eficiencia en el sitio objetivo (≥ 60%) y efectos fuera del sitio objetivo mínimos, según lo predicho por la herramienta. Los sgARN seleccionados se clonaron posteriormente en el vector pX333, que coexpresa Cas9 y los sgARN para una edición eficiente del genoma. El plásmido se transfectó en células HCT116 utilizando protocolos de transfección optimizados. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 μg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado. Luego, se recolectaron las células transfectadas y se sembraron aproximadamente 500 células en una placa de 10 cm para facilitar la expansión clonal.

Las células se cultivaron durante 10 días para permitir la formación de colonias monoclonales. Las colonias que alcanzaron un diámetro de 1 mm se aislaron manualmente y se transfirieron a placas de 48 pocillos para una mayor expansión. Una vez que alcanzaron la confluencia, se extrajo el ADN genómico de cada clon y la región objetivo se amplificó mediante PCR. El KO exitoso del fragmento que contiene el SNP se confirmó mediante secuenciación de Sanger, lo que validó la eliminación precisa de la región objetivo.

Mutación puntual

Para delimitar con precisión las consecuencias funcionales de los SNP en la actividad de los potenciadores, implementamos una estrategia de reparación dirigida por homología (HDR) basada en CRISPR-Cas9 para introducir mutaciones puntuales específicas en los sitios de los SNP. Este enfoque permitió la generación de líneas celulares isogénicas con alteraciones genéticas definidas, proporcionando una plataforma sólida para investigar el papel de los SNP en los mecanismos de regulación génica. Los sgARN dirigidos a los sitios de los SNP se diseñaron utilizando la herramienta en línea CHOPCHOP (https://chopchop.cbu.uib.no/) y se clonaron en el vector pSpCas9(BB)-2A-Puro (PX459) V2.0 (Addgene, n.º 62988). Este vector coexpresa Cas9 y el sgARN y confiere resistencia a la puromicina para una selección eficiente de las células transfectadas. Se sintetizó in vitro un oligodesoxinucleótido de cadena simple (ssODN) de 99 nt para que sirviera como plantilla de reparación para la HDR. El ssODN se diseñó para incorporar la mutación puntual deseada, al tiempo que se introducen mutaciones silenciosas en el sitio adyacente al espaciador del protómero para evitar una nueva escisión por Cas9. Este diseño aseguró una edición precisa y minimizó las alteraciones genómicas no deseadas. El plásmido pX459 (500 ng) y el donante ssODN (20 μM) se cotransfectaron en células 116 en placas de 12 pocillos utilizando los protocolos JetPRIME. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 μg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado, lo que aseguró el enriquecimiento de las células transfectadas con éxito. La selección y el clonado siguieron el protocolo de KO. Todas las secuencias de sgARN se enumeran en la tabla S16.

Metilación dirigida del ADN mediada por CRISPR-dCas9

La metilación dirigida del ADN del locus cg08640619 se logró utilizando un sistema CRISPR-dCas9-DNMT3A. En resumen, el ARN guía (sgARN) diseñado para dirigirse a la secuencia genómica que rodea a cg08640619 se clonó en un plásmido que expresa una dCas9 catalíticamente fusionada al dominio catalítico de DNMT3A. Se utilizó un sgARN no dirigido como control negativo. Las células HCT116 se recolectaron 72 horas después de la transfección. El impacto en la expresión génica se evaluó mediante análisis de transcripción inversa (RT)-qPCR de los genes diana putativos KIRREL1 y ETV3.

Ensayo CCK-8

Las células de los grupos de control y experimentales se sembraron en placas de 96 pocillos a una densidad de 1 × 104 células por pocillo, con cuatro pocillos de réplica para cada grupo. En los momentos designados (0, 24, 48 y 72 horas), se añadió el reactivo Cell Counting Kit-8 (CCK-8) [UElandy (UE), 10 μl por pocillo] a cada pocillo. Después de 1,5 horas de incubación a 37 °C, se midió la absorbancia a 450 nm utilizando un lector de microplacas (BioTek). La tasa de proliferación relativa se calculó normalizando los valores de OD450 (densidad óptica a 450 nm) de los grupos experimentales con respecto a los del grupo de control a las 12 horas (tasa de proliferación relativa = OD450 grupo experimental/OD450 grupo de control).

Ensayo de formación de colonias

Se sembraron células (500 por pocillo) en placas de seis pocillos. Después de 14 días, las colonias se fijaron [4% de paraformaldehído (PFA)], se tiñeron (0,1% de cristal violeta) y se contaron (ImageJ; colonias de > 50 células). El experimento se realizó por triplicado.

Ensayo de curación de heridas

El ensayo de curación de heridas es un método simple y ampliamente utilizado para evaluar la capacidad de migración celular mediante la creación de un arañazo artificial en un monocapa celular confluente y el seguimiento de la migración de las células para cerrar la herida. Las células se cultivaron en placas de seis pocillos hasta alcanzar una confluencia del 90%, seguido de la creación de arañazos lineales uniformes utilizando la punta de una pipeta de 10 μl. Después del lavado con solución salina tamponada con fosfato (PBS), el medio se reemplazó con medio fresco que contenía 2% de suero fetal bovino (SFB). El cierre de la herida se documentó a las 0 y 72 horas utilizando un microscopio Nikon Eclipse Ti. La tasa de migración se analizó cuantitativamente con el software ImageJ y se calculó de la siguiente manera: Tasa de migración (%) = [(Área inicial - Área final)/Área inicial] × 100%.

Ensayo Transwell

El ensayo Transwell se utiliza para investigar las capacidades de migración e invasión celular. Las células se cultivaron hasta alcanzar una confluencia del 80 al 90% y se privaron de suero durante 8 horas antes de la tripsinización. Las células desprendidas se resuspendieron en medio libre de suero y se cargaron cuidadosamente 1 × 105 células en suspensión en la cámara superior. La cámara inferior se llenó con medio que contenía 10% de SFB como agente quimiotáctico. Después de una incubación de 48 horas a 37 °C con 5% de CO2, las células no migratorias y el Matrigel residual se eliminaron de la cámara superior utilizando hisopos de algodón. La membrana se lavó dos veces con PBS frío, se fijó con 0,4% de paraformaldehído durante 10 minutos y se tiñó con 0,1% de cristal violeta durante 10 minutos. Las células invasoras se cuantificaron mediante examen microscópico y software de análisis de imágenes.

Silenciamiento génico mediante siRNA y shRNA

El silenciamiento génico se realizó utilizando oligonucleótidos siRNA dirigidos a los genes de interés. Las células se sembraron en placas de seis pocillos a 2,5 × 105 células por pocillo y se cultivaron durante 24 horas para alcanzar una confluencia del 40 al 50%. Los complejos de transfección se prepararon incubando 30 nM de siRNA (concentración final) con 4 μl de reactivo de transfección JetPrime (Polyplus) en 200 μl de medio libre de suero durante 10 minutos a temperatura ambiente. Los complejos se añadieron por goteo a las células en medio completo. Después de 48 horas de incubación, las células se recolectaron para los ensayos funcionales o la extracción de ARN.

Las secuencias de shRNA se diseñaron y se clonaron en el vector pLKO.1 (Addgene, n.º 10878). Las partículas lentivirales se empaquetaron en células HEK293T utilizando los plásmidos psPAX2 y pMD2.G. Las células se infectaron con las partículas lentivirales empaquetadas y se seleccionaron con puromicina (2 μg/ml) durante 72 horas. La eficiencia del silenciamiento (> 70% de reducción del ARNm) se validó mediante RT-qPCR. Todas las secuencias de siRNA y shRNA se enumeran en la tabla S17.

Inmunoprecipitación de cromatina

Los ensayos ChIP se llevaron a cabo como se describió anteriormente. En resumen, las células se fijaron con 1% de formaldehído (Sigma-Aldrich, F8775) durante 10 minutos a temperatura ambiente y se utilizó una concentración final de 125 mM de glicina para apagar la reacción. Las células se lavaron dos veces con PBS frío y se recolectaron en PBS frío mediante raspado y luego se resuspendieron en tampón de lisis hipótónico [20 mM de Hepes (pH 7,9) con 10 mM de KCl, 10% de glicerol, 1 mM de ditiotreitol y un cóctel de inhibidores de proteasas (04693124001, Roche)]. El pellet de núcleos se resuspendió en tampón de ensayo de inmunoprecipitación de radio (RIPA) [10 mM de tris-HCl (pH 8,0) con 140 mM de NaCl, 1% de Triton X-100, 1 mM de EDTA, 0,1% de SDS, 0,1% de desoxicolato de sodio e inhibidor de proteasas]. Los extractos nucleares se sonicaron para generar fragmentos de cromatina en un sonicador enfocado Covaris M220 (inserto, microTUBE de 130 μl; temperatura, 7 °C; potencia incidente máxima, 75 W; factor de ciclo de trabajo, 5%; tiempo de tratamiento, 300 s). La cromatina se limpió mediante centrifugación a 14 000 g durante 10 minutos a 4 °C y se diluyó con 1 ml de tampón RIPA por reacción. Luego, la cromatina se preclarificó con Sepharose 4 Fast Flow de proteína G (17061801, GE) durante 2 horas y se reservaron 30 μl de muestra como entrada y se sometieron a inmunoprecipitación con anticuerpo durante la noche. Al día siguiente, las cuentas de agarosa Pierce Protein A/G Plus (Thermo Fisher Scientific, 20423) que se bloquearon con RIPA que contenía 0,5% de albúmina sérica bovina se añadieron a las muestras, seguido de incubación durante 2 horas en una sala fría con agitación. Luego, las cuentas se lavaron tres veces con tampón RIPA, seguido de dos veces en tampón RIPA con 0,5 M de NaCl, una vez en tampón de LiCl [250 mM de LiCl, 1 mM de EDTA, 0,5% de IGEPAL CA-630, 0,1% de desoxicolato de sodio y 10 mM de tris-HCl (pH 8,0)] y dos veces en tampón TE frío (TE). Cada vez, las cuentas se agitaron suavemente durante 5 minutos. Después del lavado, tanto las cuentas como la muestra de entrada se añadieron con 150 μl de tampón de extracción (1% de SDS en 1× TE, 12 μl de 5 M de NaCl y 10 μg de RNasa A) y se incubaron a 37 °C durante 1 hora con agitación. Los fragmentos de ADN se purificaron utilizando columnas de purificación de ADN para la posterior qPCR con cebadores enumerados en la tabla S18.

La preparación de la biblioteca para los SNP heterocigotos implicó un protocolo de PCR de dos pasos. La primera PCR amplificó los fragmentos diana que flanquean el locus del SNP, incorporando adaptadores de secuenciación universales. La segunda PCR añadió un índice específico de la muestra utilizando los amplicones de la primera PCR como plantilla. La preparación de la biblioteca para el ChIP-seq de H3K27ac se realizó utilizando el kit de preparación de ADN NEBNext (NEB, E7645). La biblioteca se secuenció utilizando la plataforma DNBSEQ-T7 PE150 (Plataforma de servicios públicos biofarmacéuticos, Nanjing).

qPCR en tiempo real

Se extrajo el ARN total de las células utilizando el reactivo TRIzol y se cuantificó. El ADN genómico se eliminó tratando 1,5 μg de ARN con 4 μl de la mezcla 4× gDNA wiper en un tubo de PCR de 200 μl, ajustando el volumen final a 16 μl con agua libre de RNasa. La mezcla se mezcló a fondo por pipeteo y se incubó a 42 °C durante 2 minutos. La transcripción inversa se realizó añadiendo 4 μl de 5× HiScript II qRT SuperMix II, seguido de la mezcla e incubación a 50 °C durante 15 minutos y a 85 °C durante 5 segundos. El cDNA resultante se diluyó 20 veces con agua y se mezcló bien para su uso inmediato en qPCR. El gliceraldehído-3-fosfato deshidrogenasa sirvió como gen de referencia interno y la expresión génica relativa se calculó utilizando el método 2-ΔΔCt. Todos los cebadores de qPCR se enumeran en la tabla S17.

Métodos analíticos

Definición de potenciadores y construcción de un mapa de potenciadores a nivel del pán-cáncer

Seleccionamos líneas celulares representativas de 20 tipos de cáncer comunes. Para cada línea celular, se obtuvieron conjuntos de datos de ChIP-seq de H3K27ac y ATAC-seq o secuenciación de desoxirribonucleasa (DNase-seq) del proyecto ENCODE o del Banco de datos de expresión génica (GEO) (las fuentes de datos se enumeran en la tabla S2). Para las líneas celulares con múltiples réplicas biológicas, los picos se integraron utilizando el marco de la tasa de descubrimiento irreproducible (IDR) (^[53]). Se aplicó un umbral de IDR de 0,05 siguiendo las pautas de ENCODE (^[54]) para retener los picos de alta confianza que exhiben ≥ 95% de reproducibilidad entre las réplicas (una fracción irreproducible < 5%), lo que garantiza que solo se incluyan señales regulatorias sólidas y consistentes. Para las líneas celulares que no están incluidas en ENCODE, se descargaron los archivos de picos procesados de ChIP-Atlas, que aplica una canalización de análisis estandarizada: alineación de lecturas utilizando Bowtie2 (versión 2.4.5), eliminación de duplicados de PCR con SAMtools (versión 1.17) y llamada de picos mediante MACS2 (versión 2.2.9.1) con un umbral de significación de q < 1 × 10-5.

Por último, la intersección de los picos de H3K27ac y ATAC/DNasa se definió como las posibles regiones potenciadoras para cada línea celular, asegurando que estuvieran presentes tanto las señales de modificación de histonas como de accesibilidad de la cromatina. Para evaluar sistemáticamente la actividad potenciadora específica de cada tejido, anotamos un total de 18.880 regiones potenciadoras de CRC. Utilizando bedtools intersect (versión 2.30.0), cada región potenciadora de CRC se comparó con las regiones potenciadoras de los otros 19 tipos de cáncer, y se calculó el número de superposiciones para cada potenciador. Todos los análisis se realizaron en el genoma de referencia hg19 (GRCh37).

Análisis de datos STARR-seq

Los fragmentos de STARR-seq se alinearon a una biblioteca de oligonucleótidos personalizada utilizando Bowtie2 (versión 2.5.2). Se determinaron los recuentos de fragmentos para cada oligonucleótido por separado para las bibliotecas de ADN de entrada y ARN de salida, generando las matrices de recuento correspondientes. La normalización se realizó utilizando el paquete DESeq2 (versión 1.40.0) en R. Para cada oligonucleótido, se calculó una puntuación de actividad potenciadora (EAS) de la siguiente manera: EAS = normalizado (ARN + 1) / normalizado (ADN). Para cada SNP o sitio de metilación, se calculó el cambio en la actividad potenciadora (EASalt/EASref o EASmetilado/EASno metilado). Para definir esto geométricamente, primero ordenamos los valores de log2FC en orden ascendente. Para el umbral del lado derecho, restablecemos todos los valores menores que 0 a cero. Luego, definimos el punto de inflexión como el punto donde una línea recta con pendiente es tangente a la curva [pendiente = máx(log2FC) ? mín(log2FC)]. De manera similar, para el umbral del lado izquierdo, restablecemos todos los valores mayores que 0 a cero y utilizamos el mismo cálculo de pendiente para determinar el punto de inflexión. Los SNPs o sitios de metilación que exhiben valores de log2FC mayores que el umbral del lado derecho o menores que el umbral del lado izquierdo se seleccionan posteriormente para el análisis diferencial. La significación estadística de las diferencias en la actividad potenciadora se evaluó utilizando una prueba t de dos colas, con un valor de P < 0,05 considerado significativo.

Análisis de datos RNA-seq

Los fragmentos de RNA-seq se alinearon al genoma de referencia humano (hg19) utilizando HISAT2 (versión 2.2.1). Los recuentos de fragmentos a nivel de gen se cuantificaron utilizando featureCounts (versión 2.0.1). El análisis de expresión génica diferencial se realizó utilizando el paquete DESeq2 (versión 1.40.0). Los genes con log2FC > 0,75 y valor de P ajustado (Padj) < 0,05 se consideraron DEG. La visualización de los DEG se realizó utilizando el paquete ggplot2 (versión 3.5.1), y el análisis de enriquecimiento GO se realizó utilizando clusterProfiler (versión 4.8.3). Los conjuntos de datos de RNA-seq para las líneas celulares SW480 y SW620 se obtuvieron de GEO (www.ncbi.nlm.nih.gov/geo/query/acc.cgi acc=GSE131948).

Análisis de enriquecimiento de motivos de HOMER

Para identificar los factores de transcripción (FT) potencialmente asociados con los cambios en la actividad potenciadora, utilizamos el comando “findMotifsGenome.pl” en el software HOMER (versión 4.11). Los potenciadores diferenciales identificados a partir del ensayo STARR-seq se utilizaron como el conjunto de entrada y las regiones genómicas seleccionadas aleatoriamente sirvieron como fondo. Todos los análisis se realizaron utilizando el genoma de referencia hg19.

Escaneo de motivos con FIMO

Construimos archivos de secuencia en formato FASTA que contienen las regiones genómicas que rodean cada SNP objetivo. Se descargaron archivos de motivos de FT no redundantes para eucariotas de la base de datos JASPAR CORE (JASPAR, una base de datos de perfiles de unión de FT). El escaneo de motivos se realizó utilizando la herramienta en línea FIMO (Find Individual Motif Occurrences; https://meme-suite.org/meme/tools/fimo), que identifica las coincidencias de motivos estadísticamente significativas dentro de las secuencias de entrada.

Extracción de datos de eQTL

Para evaluar de manera integral el potencial regulador de los SNPs identificados, integramos la evidencia de eQTL tanto de contextos normales como tumorales. Contexto de tejido normal: Las estadísticas resumidas para los eQTL en el tejido normal del colon se obtuvieron a través de la plataforma FIVEx (^[55]) (https://fivex.sph.umich.edu/), que agrega y sirve datos procesados del Catálogo de eQTL del Instituto Europeo de Bioinformática (EBI) (www.ebi.ac.uk/eqtl/). Para cada variante objetivo, consultamos el punto final de la API pública de FIVEx (https://fivex.sph.umich.edu/) y recuperamos los resultados en formato JSON. Los resultados se filtraron para retener solo los registros del colon sigmoide (colonsigmoid) y el colon transverso (colontransverse).

Contexto tumoral primario: Los conjuntos de datos de eQTL de TCGA-COAD y adenocarcinoma rectal (READ) se obtuvieron de PancanQTL (^[56]) (https://gonglab.hzau.edu.cn/PancanQTL). Se consideró que un SNP tenía evidencia de eQTL de respaldo si estaba significativamente asociado con la expresión génica en el conjunto de datos de tejido normal o en los conjuntos de datos de tumores primarios (P_ < 0,05, FDR < 0,05). Este enfoque inclusivo nos permitió capturar los SNPs con un amplio potencial regulador en diferentes estados fisiológicos.

Construcción y evaluación de 14 modelos de clasificación de sondas de metilación de cáncer

Obtuvimos los datos de metilación TCGA 450k de UCSC Xena (https://xenabrowser.net/datapages/). Después del ordenamiento y el análisis estadístico, identificamos 14 tipos de cáncer con al menos 10 muestras normales y tumorales cada uno. Extraímos los datos del nivel de metilación de cg08640619 y cg25982657, junto con las etiquetas de las muestras, de los datos de metilación de estos tipos de cáncer. Utilizando la biblioteca Python scikit (versión 3.12.10), construimos conjuntos de datos de entrenamiento y prueba con la función traintestsplit, estableciendo el parámetro testsize en 0,4 para garantizar una representación equilibrada de las muestras normales y tumorales. Implementamos dos algoritmos de clasificación: LogR y SVM. Para evaluar y comparar su rendimiento durante la fase de selección del modelo, utilizamos la validación cruzada de 10 pliegues en el conjunto de entrenamiento y evaluamos la precisión de la clasificación utilizando la función crossval_score de scikit-learn. Luego, el AUROC, calculado a través del módulo de métricas de scikit-learn, se utilizó como la métrica principal para evaluar el rendimiento de generalización de los clasificadores seleccionados en datos no vistos. Por último, resumimos y representamos los valores de AUROC de los clasificadores LogR y SVM construidos con cg08640619 y cg25982657 para los 14 tipos de cáncer en gráficos de barras.

Análisis estadístico experimental

El análisis estadístico se realizó utilizando GraphPad Prism (versión 8, GraphPad Software Inc.). Los experimentos se repitieron al menos tres veces, y los resultados se presentan como medias ± DE. Se utilizaron pruebas t o análisis de varianza (ANOVA) para las comparaciones estadísticas, con valores de P < 0,05 considerados estadísticamente significativos.

Diseño de oligonucleótidos

Establecimos un conjunto integral de SNPs a través de criterios integrados de relevancia funcional y de enfermedad: (i) aplicamos un filtrado estricto de MAF [0,05 en AMR/EUR/EAS/SAS en el Proyecto de los 1000 Genomas (^[15])]; (ii) ampliamos 238 variantes de GWAS de CRC mediante la imputación de proxy de LD (R2 ≥ 0,5, ±500 kb), lo que resultó en 6125 SNPs adicionales; y (iii) priorizamos de forma independiente 24.526 SNPs localizados en potenciadores definidos por la coaccesibilidad de la cromatina [intersección de H3K27ac ChIP-seq (ENCSR661KMA/ENCSR000EUT) y picos de ATAC-seq (ENCSR872WGW) en células HCT-116]. El conjunto de unión comprendió 30.790 SNPs únicos (6363 derivados de GWAS y 24.526 basados en potenciadores, con 99 superposiciones).

Para los oligonucleótidos específicos de metilación, seleccionamos (i) los CpGs del núcleo del potenciador: más de 134.000 dinucleótidos CpG dentro de ±60 pb de los picos de H3K27ac; y (ii) los CpGs diferenciales de CRC (| cambio | > 0,3, FDR < 0,05): 908 sondas localizadas en potenciadores de las matrices Illumina 450K de la cohorte TCGA-CRC, GSE139404, GSE77955, GSE53051, GSE48684, GSE42752 y GSE40055. El diseño de los oligonucleótidos se adaptó a un sistema de secuenciación de Beijing Genomics Institute (BGI), y las dos bibliotecas se sintetizaron por separado por CustomArray. En resumen, cada oligonucleótido contiene 120 pb de secuencia genómica que encierra los sitios de SNP o CpG y secuencias de flanqueo constantes (ascendente, 5'-ATGGCTACGATCCGACTT; y descendente, AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA-3') en ambos extremos, que se utilizaron para el clonado.

Preparación del plásmido SNP-STARR-seq

El vector hORI-STARR-seq fue proporcionado por A. Stark (Instituto de Investigación de Patología Molecular, Viena, Austria). La biblioteca de SNPs amplificada por PCR (20 ng) se insertó en los sitios Age I/Sal I del plásmido (100 ng) utilizando un clonado basado en recombinación altamente eficiente (ClonExpress, Vazyme, China). Para evitar sesgos durante el clonado, realizamos un total de 20 reacciones de recombinación separadas y agrupamos cada cuatro reacciones. Después de 1 × purificación de ADN con cuentas AMPure XP, utilizamos los vectores que contienen SNPs (2 μl) para transformar Escherichia coli DH5α (DH5α; 50 μl). Nuevamente, realizamos 20 reacciones de transformación separadas y agrupamos cada cuatro transformaciones. Luego, los plásmidos se extrajeron y purificaron a una alta concentración (>1 μg/μl). Se realizó una amplificación por PCR a partir de los plásmidos y se secuenció como control de entrada en ciclos de 2 × 100 de extremos emparejados.

Preparación del plásmido Metil-STARR-seq

pmSTARRseq1 fue un regalo de J. Tung. Los métodos de clonado son similares a los del SNP-STARR-seq mencionado anteriormente, excepto que la biblioteca de metilación amplificada por PCR (10 ng) se insertó en los sitios Spe I/Nco I del plásmido (100 ng). Los plásmidos se transformaron en células GT115 extra-competentes (con selección con zeocina a 50 μg/ml). Luego, metilamos la mitad de la biblioteca de plásmidos con M.SssI [ADN, 4 U/μg; New England Biolabs (NEB)] y tratamos la otra mitad como control simulado (utilizando el mismo protocolo de reacción, pero reemplazando la enzima con agua destilada). Para confirmar el estado de metilación del ADN de los plásmidos, primero digerimos el ADN del plásmido (1 μg) incubándolo (a 37 °C) con Xba I y Nhe I en un total de cinco reacciones. Se añadió aproximadamente 1 ng de fragmento de proteína fluorescente verde no metilada a cada muestra para evaluar la eficiencia de la conversión de bisulfito. Realizamos reparación de extremos, adición de cola A y ligadura de adaptadores. Luego, las muestras se sometieron a conversión de bisulfito de sodio utilizando un kit EZ DNA Methylation-Gold (Zymo Research). Las bibliotecas se amplificaron por PCR (KAPA HiFi Uracil+ DNA Polymerase) y se secuenciaron en la plataforma Illumina.

Cultivo celular y transfección

Cada grupo de oligonucleótidos se transfectó en células embrionarias humanas (HEK 293T) o líneas celulares de CRC (HCT116, SW480 y SW620) utilizando jetOPTIMUS (Polyplus). Las células se cultivaron en condiciones normales con 5 % de CO2 a 37 °C. Específicamente, se mezclaron 2 μg de plásmidos con 3 μl de reactivos de transfección para la transfección en células cultivadas en un pozo de una placa de seis pocillos. Agrupamos cada seis pocillos (~2 × 107 células) para una repetición biológica.

Extracción de ARNm y secuenciación

Treinta y seis horas después de la transfección, se extrajo el ARN total con un kit RNeasy (TAKARA) y se enriqueció el ARNm con el módulo de aislamiento de ARNm poli(A) magnético (NEB, E7490L). La primera cadena de ADN complementario se sintetizó utilizando un cebador específico (5'-CTCATCAATGTATCTTATCATGTCTG para SNP-STARR-seq y 5'-CAAACTCATCAATGTATCTTATCATG para metil-STARR-seq) con SMARTScribe Reverse Transcriptase (Takara Bio, 639538). Amplificamos el ADN complementario obtenido de la transcripción inversa para la secuenciación de BGI utilizando una PCR anidada de dos pasos. Para la PCR de unión (15 ciclos), el cebador ascendente abarca la unión del intrón sintético y amplifica específicamente el ADN complementario del reportero sin amplificar ningún ADN plasmídico (ver cebadores en la tabla S12). El producto purificado total sirvió como plantilla para la PCR de índice (10 ciclos). El ADN se purificó con cuentas AMPure y se secuenció en ciclos de 2 × 100 de extremos emparejados con el secuenciador BGI. En total, se realizaron tres repeticiones biológicas con tres repeticiones técnicas cada una para todas las líneas celulares.

Ensayo de reportero de luciferasa dual

SPANISH TRANSLATION:

Para investigar el impacto funcional de los SNPs o los sitios de metilación en la actividad de los potenciadores, se realizaron ensayos de reportero de luciferasa dual. En resumen, fragmentos de ADN de 120 pares de bases que contienen distintos alelos de SNPs o fragmentos que contienen CpG se clonaron en el vector reportero de luciferasa pGL4.23 (sitios Kpn I/Xho I) o en el vector pCpG-free-lucia (sitios BamH I/Spe I), respectivamente. El plásmido recombinante pGL4.23 o pCpG-free-lucia (reportero de luciferasa de luciérnaga) se cotransfectó con el vector pGL4.74 para normalizar las variaciones en la eficiencia de la transfección y la viabilidad celular. Después de 48 horas de transfección, se midió la actividad de la luciferasa utilizando el kit de ensayo de luciferasa Dual-Lumi (Beyotime, RG088S). Actividad relativa del potenciador = Luminiscencia de luciérnaga/Luminiscencia de renilla. Los cebadores de clonación, así como las secuencias de los sitios variantes detectados, se enumeran en la tabla S13. Los datos originales de la luminiscencia de luciérnaga y renilla se enumeran en la tabla S14.

Secuenciación de ARN

Se extrajo el ARN total de las células utilizando el kit de extracción de ARN universal MiniBEST de TAKARA. Se enriqueció el ARN poliadenilado a partir de 5 µg de ARN total, seguido de una elución directa utilizando un sistema de transcripción inversa preformulado para la síntesis de ADNc. El ARN residual se eliminó mediante tratamiento con ribonucleasa (RNasa) A/H. El ADNc se purificó y recuperó utilizando cuentas AMPure XP, y se ligaron adaptadores en el extremo 5'. Después de dos rondas de purificación con cuentas AMPure XP, el ADN se amplificó utilizando los cebadores de secuenciación BGI (BGI-system F y BGI-Barcode R), con distintos cebadores posteriores utilizados para diferentes muestras para permitir el multiplexado de muestras. Tras la amplificación, el ADN se enriqueció aún más utilizando cuentas AMPure XP antes de agruparlo para la secuenciación.

Genotipado

Se extrajo el ADN genómico de las células siguiendo el protocolo del fabricante del kit de extracción de ADN. Se diseñaron pares de cebadores que flanquean el locus de interés del SNP (tabla S15) y se utilizaron para la amplificación por PCR con la polimerasa Taq. Los productos de PCR se enviaron directamente a una empresa de biotecnología para la secuenciación de Sanger. El genotipo del locus de SNP se determinó analizando los picos del cromatograma de secuenciación en la posición objetivo (SnapGene).

Edición génica mediada por CRISPR-Cas9

Inactivación

Para dilucidar el papel funcional de regiones específicas que contienen SNPs como posibles potenciadores, utilizamos la edición génica mediada por CRISPR-Cas9 para generar la inactivación dirigida de los fragmentos del sitio de SNP. Este enfoque permitió la evaluación de la actividad del potenciador mediante el análisis cuantitativo de los cambios en la expresión génica posterior.

Se diseñaron dos ARN guía simples (sgARN) que flanquean la región del SNP objetivo utilizando la plataforma de diseño de CRISPR CHOPCHOP (https://chopchop.cbu.uib.no/). Los sgARN se seleccionaron en base a criterios estrictos, que incluyen una alta eficiencia en el sitio objetivo (≥ 60%) y efectos fuera del sitio mínimos, según lo predicho por la herramienta. Los sgARN seleccionados se clonaron posteriormente en el vector pX333, que coexpresa Cas9 y los sgARN para una edición eficiente del genoma. El plásmido se transfectó en células HCT116 utilizando protocolos de transfección optimizados. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 µg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado. Las células transfectadas se cosecharon y se sembraron aproximadamente 500 células en una placa de 10 cm para facilitar la expansión clonal.

Las células se cultivaron durante 10 días para permitir la formación de colonias monoclonales. Las colonias que alcanzaron un diámetro de 1 mm se aislaron manualmente y se transfirieron a placas de 48 pocillos para una mayor expansión. Una vez que se alcanzó la confluencia, se extrajo el ADN genómico de cada clon y la región objetivo se amplificó por PCR. La inactivación exitosa del fragmento que contiene el SNP se confirmó mediante la secuenciación de Sanger, lo que validó la eliminación precisa de la región objetivo.

Mutación puntual

Para delinear con precisión las consecuencias funcionales de los SNPs en la actividad del potenciador, implementamos una estrategia de reparación dirigida por homología (HDR) basada en CRISPR-Cas9 para introducir mutaciones puntuales específicas en los sitios de SNP. Este enfoque permitió la generación de líneas celulares isogénicas con alteraciones genéticas definidas, proporcionando una plataforma sólida para investigar el papel de los SNPs en los mecanismos de regulación génica. Los sgARN dirigidos a los sitios de SNP se diseñaron utilizando la herramienta en línea CHOPCHOP (https://chopchop.cbu.uib.no/) y se clonaron en el vector pSpCas9(BB)-2A-Puro (PX459) V2.0 (Addgene, n.º 62988). Este vector coexpresa Cas9 y el sgARN y confiere resistencia a la puromicina para una selección eficiente de las células transfectadas. Se sintetizó in vitro un oligonucleótido de cadena simple (ssODN) de 99 nt para servir como plantilla de reparación para HDR. El ssODN se diseñó para incorporar la mutación puntual deseada, al tiempo que se introducen mutaciones silenciosas en el sitio adyacente al espaciador del protómero para evitar una nueva escisión por Cas9. Este diseño aseguró una edición precisa y minimizó las alteraciones genómicas no deseadas. El plásmido pX459 (500 ng) y el donante de ssODN (20 µM) se cotransfectaron en células 116 en placas de 12 pocillos utilizando los protocolos JetPRIME. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 µg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado, lo que aseguró el enriquecimiento de las células transfectadas con éxito. La selección y el clonado siguieron el protocolo de inactivación. Todas las secuencias de sgARN se enumeran en la tabla S16.

Inactivación

Para dilucidar el papel funcional de regiones específicas que contienen SNPs como posibles potenciadores, utilizamos la edición génica mediada por CRISPR-Cas9 para generar la inactivación dirigida de los fragmentos del sitio de SNP. Este enfoque permitió la evaluación de la actividad del potenciador mediante el análisis cuantitativo de los cambios en la expresión génica posterior.

Se diseñaron dos ARN guía simples (sgARN) que flanquean la región del SNP objetivo utilizando la plataforma de diseño de CRISPR CHOPCHOP (https://chopchop.cbu.uib.no/). Los sgARN se seleccionaron en base a criterios estrictos, que incluyen una alta eficiencia en el sitio objetivo (≥ 60%) y efectos fuera del sitio mínimos, según lo predicho por la herramienta. Los sgARN seleccionados se clonaron posteriormente en el vector pX333, que coexpresa Cas9 y los sgARN para una edición eficiente del genoma. El plásmido se transfectó en células HCT116 utilizando protocolos de transfección optimizados. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 µg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado. Las células transfectadas se cosecharon y se sembraron aproximadamente 500 células en una placa de 10 cm para facilitar la expansión clonal.

Las células se cultivaron durante 10 días para permitir la formación de colonias monoclonales. Las colonias que alcanzaron un diámetro de 1 mm se aislaron manualmente y se transfirieron a placas de 48 pocillos para una mayor expansión. Una vez que se alcanzó la confluencia, se extrajo el ADN genómico de cada clon y la región objetivo se amplificó por PCR. La inactivación exitosa del fragmento que contiene el SNP se confirmó mediante la secuenciación de Sanger, lo que validó la eliminación precisa de la región objetivo.

Mutación puntual

Para delinear con precisión las consecuencias funcionales de los SNPs en la actividad del potenciador, implementamos una estrategia de reparación dirigida por homología (HDR) basada en CRISPR-Cas9 para introducir mutaciones puntuales específicas en los sitios de SNP. Este enfoque permitió la generación de líneas celulares isogénicas con alteraciones genéticas definidas, proporcionando una plataforma sólida para investigar el papel de los SNPs en los mecanismos de regulación génica. Los sgARN dirigidos a los sitios de SNP se diseñaron utilizando la herramienta en línea CHOPCHOP (https://chopchop.cbu.uib.no/) y se clonaron en el vector pSpCas9(BB)-2A-Puro (PX459) V2.0 (Addgene, n.º 62988). Este vector coexpresa Cas9 y el sgARN y confiere resistencia a la puromicina para una selección eficiente de las células transfectadas. Se sintetizó in vitro un oligonucleótido de cadena simple (ssODN) de 99 nt para servir como plantilla de reparación para HDR. El ssODN se diseñó para incorporar la mutación puntual deseada, al tiempo que se introducen mutaciones silenciosas en el sitio adyacente al espaciador del protómero para evitar una nueva escisión por Cas9. Este diseño aseguró una edición precisa y minimizó las alteraciones genómicas no deseadas. El plásmido pX459 (500 ng) y el donante de ssODN (20 µM) se cotransfectaron en células 116 en placas de 12 pocillos utilizando los protocolos JetPRIME. Después de 24 horas, se reemplazó el medio y las células se sometieron a selección con puromicina (1,5 µg/ml). La selección se mantuvo hasta que se observó una muerte celular casi completa en el grupo de control no transfectado, lo que aseguró el enriquecimiento de las células transfectadas con éxito. La selección y el clonado siguieron el protocolo de inactivación. Todas las secuencias de sgARN se enumeran en la tabla S16.

Edición dirigida de la metilación del ADN mediada por CRISPR-dCas9

La metilación dirigida del ADN del locus cg08640619 se logró utilizando un sistema CRISPR-dCas9-DNMT3A. En resumen, se clonó un ARN guía (sgARN) diseñado para dirigirse a la secuencia genómica que rodea a cg08640619 en un plásmido que expresa una dCas9 catalíticamente fusionada al dominio catalítico de DNMT3A. Se utilizó un sgARN no dirigido como control negativo. Las células HCT116 se cosecharon 72 horas después de la transfección. El impacto en la expresión génica se evaluó mediante el análisis de transcripción inversa (RT)-qPCR de los genes diana putativos KIRREL1 y ETV3.

Ensayo CCK-8

Las células de los grupos de control y experimentales se sembraron en placas de 96 pocillos a una densidad de 1 × 104 células por pocillo, con cuatro pocillos de réplica para cada grupo. En los momentos designados (0, 24, 48 y 72 horas), se añadió el reactivo Cell Counting Kit-8 (CCK-8) [UElandy (UE), 10 µl por pocillo] a cada pocillo. Después de 1,5 horas de incubación a 37 °C, se midió la absorbancia a 450 nm utilizando un lector de microplacas (BioTek). La tasa de proliferación relativa se calculó normalizando los valores de OD450 (densidad óptica a 450 nm) de los grupos experimentales con respecto a los del grupo de control a las 12 horas (Tasa de proliferación relativa = OD450 grupo experimental/OD450 grupo de control).

Ensayo de formación de colonias

Se sembraron células (500 por pocillo) en placas de seis pocillos. Después de 14 días, las colonias se fijaron [4% de paraformaldehído (PFA)], se tiñeron (0,1% de cristal violeta) y se contaron (ImageJ; colonias de > 50 células). El experimento se realizó por triplicado.

Ensayo de curación de heridas

El ensayo de curación de heridas es un método sencillo y ampliamente utilizado para evaluar la capacidad de migración celular mediante la creación de un arañazo artificial en un monocapa celular confluente y el seguimiento de la migración de las células para cerrar la herida. Las células se cultivaron en placas de seis pocillos hasta alcanzar una confluencia del 90%, seguida de la creación de arañazos lineales uniformes utilizando una punta de pipeta de 10 µl. Después del lavado con solución salina tamponada con fosfato (PBS), el medio se reemplazó con medio fresco que contenía un 2% de suero fetal bovino (SFB). El cierre de la herida se documentó a las 0 y 72 horas utilizando un microscopio Nikon Eclipse Ti. La tasa de migración se analizó cuantitativamente con el software ImageJ y se calculó como: Tasa de migración (%) = [(Área inicial - Área final)/Área inicial] × 100%.

Ensayo Transwell

El ensayo Transwell se utiliza para investigar las capacidades de migración e invasión celular. Las células se cultivaron hasta alcanzar una confluencia del 80 al 90% y se privaron de suero durante 8 horas antes de la tripsinización. Las células desprendidas se resuspendieron en un medio libre de suero y se cargaron cuidadosamente 1 × 105 células en suspensión en la cámara superior. La cámara inferior se llenó con medio que contenía un 10% de SFB como atrayente. Después de 48 horas de incubación a 37 °C con 5% de CO2, las células no migrantes y el gel de Matrigel residual se eliminaron de la cámara superior utilizando hisopos de algodón. La membrana se lavó dos veces con PBS frío, se fijó con 0,4% de paraformaldehído durante 10 minutos y se tiñó con 0,1% de cristal violeta durante 10 minutos. Las células invasoras se cuantificaron mediante el examen microscópico y el software de análisis de imágenes.

Supresión de la expresión génica mediante siRNA y shRNA

Se realizó el silenciamiento génico utilizando oligonucleótidos de siRNA dirigidos a los genes de interés. Las células se sembraron en placas de seis pocillos a una densidad de 2,5 × 105 células por pocillo y se cultivaron durante 24 horas para alcanzar una confluencia del 40 al 50 %. Los complejos de transfección se prepararon incubando 30 nM de siRNA (concentración final) con 4 μl de reactivo de transfección JetPrime (Polyplus) en 200 μl de medio libre de suero durante 10 minutos a temperatura ambiente. Los complejos se añadieron gota a gota a las células en medio completo. Después de 48 horas de incubación, las células se recogieron para los ensayos funcionales o la extracción de ARN.

Las secuencias de shRNA se diseñaron y clonaron en el vector pLKO.1 (Addgene, n.º 10878). Las partículas lentivirales se empaquetaron en células HEK293T utilizando los plásmidos psPAX2 y pMD2.G. Las células se infectaron con los lentivirus empaquetados y se seleccionaron con puromicina (2 μg/ml) durante 72 horas. La eficiencia del silenciamiento (>70 % de reducción del ARNm) se validó mediante RT-qPCR. Todas las secuencias de siRNA y shRNA se enumeran en la tabla S17.

Inmunoprecipitación de cromatina

Los ensayos de ChIP se llevaron a cabo según lo descrito previamente. Brevemente, las células se fijaron con un 1 % de formaldehído (Sigma-Aldrich, F8775) durante 10 minutos a temperatura ambiente, y se utilizó una concentración final de 125 mM de glicina para detener la reacción. Las células se lavaron dos veces con PBS frío y se recogieron en PBS frío por raspado y luego se resuspendieron en un tampón de lisis hipótónico [20 mM de Hepes (pH 7,9) con 10 mM de KCl, 10 % de glicerol, 1 mM de ditiotreitol y un cóctel de inhibidores de proteasas (04693124001, Roche)]. El precipitado de núcleos se resuspendió en un tampón de ensayo de radioinmunoprecipitación (RIPA) [10 mM de tris-HCl (pH 8,0) con 140 mM de NaCl, 1 % de Triton X-100, 1 mM de EDTA, 0,1 % de SDS, 0,1 % de desoxicolato de sodio e inhibidor de proteasas]. Los extractos nucleares se sonicaron para generar fragmentos de cromatina en un sonicador Covaris M220 Focused-ultrasonicator (inserto, microTUBE 130 μl; temperatura, 7 °C; potencia incidente máxima, 75 W; factor de ciclo de trabajo, 5 %; tiempo de tratamiento, 300 s). La cromatina se purificó mediante centrifugación a 14 000 g durante 10 minutos a 4 °C y se diluyó con 1 ml de tampón RIPA por reacción. A continuación, la cromatina se preclarificó con Protein G Sepharose 4 Fast Flow (17061801, GE) durante 2 horas, y se reservaron 30 μl de muestra como control y se sometieron a inmunoprecipitación con anticuerpo durante la noche. Al día siguiente, se añadieron las cuentas de Pierce Protein A/G Plus Agarose (Thermo Fisher Scientific, 20423) que se habían bloqueado con RIPA que contenía un 0,5 % de albúmina sérica bovina a las muestras, seguido de incubación durante 2 horas en una cámara fría con agitación. A continuación, las cuentas se lavaron tres veces con tampón RIPA, seguido de dos veces en tampón RIPA con 0,5 M de NaCl, una vez en tampón de LiCl [250 mM de LiCl, 1 mM de EDTA, 0,5 % de IGEPAL CA-630, 0,1 % de desoxicolato de sodio y 10 mM de tris-HCl (pH 8,0)] y dos veces en tampón de tris-EDTA (TE) frío. Cada vez, las cuentas se mantuvieron durante 5 minutos con una suave agitación. Después del lavado, tanto las cuentas como la muestra de control se añadieron a 150 μl de tampón de extracción (1 % de SDS en 1× TE, 12 μl de 5 M de NaCl y 10 μg de RNasa A) y se incubaron a 37 °C durante 1 hora con agitación. Los fragmentos de ADN se purificaron utilizando columnas de purificación de ADN para la posterior qPCR con los cebadores que se enumeran en la tabla S18.

La preparación de la biblioteca para los polimorfismos de un solo nucleótido (SNP) heterocigotos implicó un protocolo de PCR de dos pasos. La primera PCR amplificó los fragmentos diana que flanqueaban el locus del SNP, incorporando adaptadores de secuenciación universales. La segunda PCR añadió un índice específico de la muestra utilizando los amplicones de la primera PCR como plantilla. La preparación de la biblioteca para la ChIP-seq de H3K27ac se realizó utilizando el kit NEBNext DNA Prep Kit (NEB, E7645). La biblioteca se secuenció utilizando la plataforma DNBSEQ-T7 PE150 (Plataforma de Servicios Públicos Biofarmacéuticos, Nanjing).

qPCR en tiempo real

Se extrajo el ARN total de las células utilizando el reactivo TRIzol y se cuantificó. El ADN genómico se eliminó tratando 1,5 μg de ARN con 4 μl de 4× gDNA wiper Mix en un tubo de PCR de 200 μl, ajustando el volumen final a 16 μl con agua libre de RNasa. La mezcla se mezcló a fondo por pipeteo y se incubó a 42 °C durante 2 minutos. La transcripción inversa se realizó añadiendo 4 μl de 5× HiScript II qRT SuperMix II, seguido de la mezcla y la incubación a 50 °C durante 15 minutos y a 85 °C durante 5 segundos. El cDNA resultante se diluyó 20 veces con agua y se mezcló bien para su uso inmediato en qPCR. La gliceraldehído-3-fosfato deshidrogenasa sirvió como gen de referencia interno, y la expresión génica relativa se calculó utilizando el método 2-ΔΔCt. Todos los cebadores de qPCR se enumeran en la tabla S17.

Ensayo CCK-8

Las células de los grupos de control y experimentales se sembraron en placas de 96 pocillos a una densidad de 1 × 104 células por pocillo, con cuatro pocillos de réplica para cada grupo. En los momentos designados (0, 24, 48 y 72 horas), se añadió el reactivo Cell Counting Kit-8 (CCK-8) [UElandy (UE), 10 μl por pocillo] a cada pocillo. Después de 1,5 horas de incubación a 37 °C, se midió la absorbancia a 450 nm utilizando un lector de microplacas (BioTek). La tasa de proliferación relativa se calculó normalizando los valores de OD450 (densidad óptica a 450 nm) de los grupos experimentales con respecto a los del grupo de control a las 12 horas (tasa de proliferación relativa = OD450 del grupo experimental/OD450 del grupo de control).

Ensayo de formación de colonias

Las células (500 por pocillo) se sembraron en placas de seis pocillos. Después de 14 días, las colonias se fijaron [4 % de paraformaldehído (PFA)], se tiñeron (0,1 % de cristal violeta) y se contaron (ImageJ; colonias de >50 células). El experimento se realizó por triplicado.

Ensayo de curación de heridas

El ensayo de curación de heridas es un método sencillo y ampliamente utilizado para evaluar la capacidad de migración celular mediante la creación de un arañazo artificial en un monocapa celular y el seguimiento de la migración de las células para cerrar la herida. Las células se cultivaron en placas de seis pocillos hasta alcanzar una confluencia del 90 %, seguido de la creación de arañazos lineales uniformes utilizando la punta de una pipeta de 10 μl. Después del lavado con solución salina tamponada con fosfato (PBS), el medio se reemplazó por medio fresco que contenía un 2 % de suero fetal bovino (SFB). El cierre de la herida se documentó a los 0 y 72 horas utilizando un microscopio Nikon Eclipse Ti. La tasa de migración se analizó cuantitativamente con el software ImageJ y se calculó de la siguiente manera: Tasa de migración (%) = [(Área inicial - Área final)/Área inicial] × 100 %.

Ensayo Transwell

El ensayo Transwell se utiliza para investigar las capacidades de migración e invasión celular. Las células se cultivaron hasta alcanzar una confluencia del 80 al 90 % y se privaron de suero durante 8 horas antes de la tripsinización. Las células desprendidas se resuspendieron en un medio libre de suero, y se cargaron cuidadosamente 1 × 105 células en suspensión en la cámara superior. La cámara inferior se llenó con medio que contenía un 10 % de SFB como agente quimiotáctico. Después de 48 horas de incubación a 37 °C con 5 % de CO2, las células no migrantes y el gel de Matrigel residual se eliminaron de la cámara superior utilizando bastoncillos de algodón. La membrana se lavó dos veces con PBS frío, se fijó con un 0,4 % de paraformaldehído durante 10 minutos y se tiñó con un 0,1 % de cristal violeta durante 10 minutos. Las células invasoras se cuantificaron mediante el examen microscópico y el software de análisis de imágenes.

Silenciamiento génico mediante siRNA y shRNA

Se realizó el silenciamiento génico utilizando oligonucleótidos de siRNA dirigidos a los genes de interés. Las células se sembraron en placas de seis pocillos a una densidad de 2,5 × 105 células por pocillo y se cultivaron durante 24 horas para alcanzar una confluencia del 40 al 50 %. Los complejos de transfección se prepararon incubando 30 nM de siRNA (concentración final) con 4 μl de reactivo de transfección JetPrime (Polyplus) en 200 μl de medio libre de suero durante 10 minutos a temperatura ambiente. Los complejos se añadieron gota a gota a las células en medio completo. Después de 48 horas de incubación, las células se recogieron para los ensayos funcionales o la extracción de ARN.

Las secuencias de shRNA se diseñaron y clonaron en el vector pLKO.1 (Addgene, n.º 10878). Las partículas lentivirales se empaquetaron en células HEK293T utilizando los plásmidos psPAX2 y pMD2.G. Las células se infectaron con los lentivirus empaquetados y se seleccionaron con puromicina (2 μg/ml) durante 72 horas. La eficiencia del silenciamiento (>70 % de reducción del ARNm) se validó mediante RT-qPCR. Todas las secuencias de siRNA y shRNA se enumeran en la tabla S17.

Inmunoprecipitación de cromatina

Los ensayos de ChIP se llevaron a cabo según lo descrito previamente. Brevemente, las células se fijaron con un 1 % de formaldehído (Sigma-Aldrich, F8775) durante 10 minutos a temperatura ambiente, y se utilizó una concentración final de 125 mM de glicina para detener la reacción. Las células se lavaron dos veces con PBS frío y se recogieron en PBS frío por raspado y luego se resuspendieron en un tampón de lisis hipótónico [20 mM de Hepes (pH 7,9) con 10 mM de KCl, 10 % de glicerol, 1 mM de ditiotreitol y un cóctel de inhibidores de proteasas (04693124001, Roche)]. El precipitado de núcleos se resuspendió en un tampón de ensayo de radioinmunoprecipitación (RIPA) [10 mM de tris-HCl (pH 8,0) con 140 mM de NaCl, 1 % de Triton X-100, 1 mM de EDTA, 0,1 % de SDS, 0,1 % de desoxicolato de sodio e inhibidor de proteasas]. Los extractos nucleares se sonicaron para generar fragmentos de cromatina en un sonicador Covaris M220 Focused-ultrasonicator (inserto, microTUBE 130 μl; temperatura, 7 °C; potencia incidente máxima, 75 W; factor de ciclo de trabajo, 5 %; tiempo de tratamiento, 300 s). La cromatina se purificó mediante centrifugación a 14 000 g durante 10 minutos a 4 °C y se diluyó con 1 ml de tampón RIPA por reacción. A continuación, la cromatina se preclarificó con Protein G Sepharose 4 Fast Flow (17061801, GE) durante 2 horas, y se reservaron 30 μl de muestra como control y se sometieron a inmunoprecipitación con anticuerpo durante la noche. Al día siguiente, se añadieron las cuentas de Pierce Protein A/G Plus Agarose (Thermo Fisher Scientific, 20423) que se habían bloqueado con RIPA que contenía un 0,5 % de albúmina sérica bovina a las muestras, seguido de incubación durante 2 horas en una cámara fría con agitación. A continuación, las cuentas se lavaron tres veces con tampón RIPA, seguido de dos veces en tampón RIPA con 0,5 M de NaCl, una vez en tampón de LiCl [250 mM de LiCl, 1 mM de EDTA, 0,5 % de IGEPAL CA-630, 0,1 % de desoxicolato de sodio y 10 mM de tris-HCl (pH 8,0)] y dos veces en tampón de tris-EDTA (TE) frío. Cada vez, las cuentas se mantuvieron durante 5 minutos con una suave agitación. Después del lavado, tanto las cuentas como la muestra de control se añadieron a 150 μl de tampón de extracción (1 % de SDS en 1× TE, 12 μl de 5 M de NaCl y 10 μg de RNasa A) y se incubaron a 37 °C durante 1 hora con agitación. Los fragmentos de ADN se purificaron utilizando columnas de purificación de ADN para la posterior qPCR con los cebadores que se enumeran en la tabla S18.

La preparación de la biblioteca para los polimorfismos de un solo nucleótido (SNP) heterocigotos implicó un protocolo de PCR de dos pasos. La primera PCR amplificó los fragmentos diana que flanqueaban el locus del SNP, incorporando adaptadores de secuenciación universales. La segunda PCR añadió un índice específico de la muestra utilizando los amplicones de la primera PCR como plantilla. La preparación de la biblioteca para la ChIP-seq de H3K27ac se realizó utilizando el kit NEBNext DNA Prep Kit (NEB, E7645). La biblioteca se secuenció utilizando la plataforma DNBSEQ-T7 PE150 (Plataforma de Servicios Públicos Biofarmacéuticos, Nanjing).

qPCR en tiempo real

Se extrajo el ARN total de las células utilizando el reactivo TRIzol y se cuantificó. El ADN genómico se eliminó tratando 1,5 μg de ARN con 4 μl de 4× gDNA wiper Mix en un tubo de PCR de 200 μl, ajustando el volumen final a 16 μl con agua libre de RNasa. La mezcla se mezcló a fondo por pipeteo y se incubó a 42 °C durante 2 minutos. La transcripción inversa se realizó añadiendo 4 μl de 5× HiScript II qRT SuperMix II, seguido de la mezcla y la incubación a 50 °C durante 15 minutos y a 85 °C durante 5 segundos. El cDNA resultante se diluyó 20 veces con agua y se mezcló bien para su uso inmediato en qPCR. La gliceraldehído-3-fosfato deshidrogenasa sirvió como gen de referencia interno, y la expresión génica relativa se calculó utilizando el método 2-ΔΔCt. Todos los cebadores de qPCR se enumeran en la tabla S17.

Métodos analíticos

Definición de potenciadores y construcción de un mapa de potenciadores pancancerosos

Seleccionamos líneas celulares representativas de 20 tipos de cáncer comunes. Para cada línea celular, se obtuvieron conjuntos de datos de ChIP-seq de H3K27ac y ATAC-seq o secuenciación de desoxirribonucleasa (DNase-seq) del proyecto ENCODE o del Gene Expression Omnibus (GEO) (las fuentes de datos se enumeran en la tabla S2). Para las líneas celulares con múltiples réplicas biológicas, los picos se integraron utilizando el marco Irreproducible Discovery Rate (IDR) (^[53]). Se aplicó un umbral de IDR de 0,05 siguiendo las pautas de ENCODE (^[54]) para retener picos de alta confianza que exhiben ≥95% de reproducibilidad entre las réplicas (una fracción irreproducible 0,75 y un valor de P ajustado (Padj) < 0,05 se consideraron DEG. La visualización de DEG se realizó utilizando el paquete ggplot2 (versión 3.5.1), y el análisis de enriquecimiento de GO se realizó utilizando clusterProfiler (versión 4.8.3). Los conjuntos de datos de RNA-seq para las líneas celulares SW480 y SW620 se obtuvieron de GEO (www.ncbi.nlm.nih.gov/geo/query/acc.cgi acc=GSE131948).

Análisis de enriquecimiento de motivos de HOMER

Para identificar los factores de transcripción (FT) potencialmente asociados con los cambios en la actividad potenciadora, utilizamos el comando “findMotifsGenome.pl” en el software HOMER (versión 4.11). Los potenciadores diferenciales identificados a partir del ensayo STARR-seq se utilizaron como el conjunto de entrada y las regiones genómicas seleccionadas aleatoriamente sirvieron como fondo. Todos los análisis se realizaron utilizando el genoma de referencia hg19.

Escaneo de motivos con FIMO

Construimos archivos de secuencia en formato FASTA que contienen las regiones genómicas que rodean cada SNP objetivo. Los archivos de motivos de FT no redundantes para eucariotas se descargaron de la base de datos JASPAR CORE (JASPAR, una base de datos de perfiles de unión de FT). El escaneo de motivos se realizó utilizando la herramienta en línea FIMO (Find Individual Motif Occurrences; https://meme-suite.org/meme/tools/fimo), que identifica coincidencias de motivos estadísticamente significativas dentro de las secuencias de entrada.

Extracción de datos de eQTL

Para evaluar de manera integral el potencial regulador de los SNPs identificados, integramos la evidencia de eQTL tanto de contextos normales como tumorales. Contexto de tejido normal: Las estadísticas resumidas para los eQTL en el tejido normal del colon se obtuvieron a través de la plataforma FIVEx (^[55]) (https://fivex.sph.umich.edu/), que agrega y sirve datos procesados del Catálogo de eQTL del European Bioinformatics Institute (EBI) (www.ebi.ac.uk/eqtl/). Para cada variante objetivo, consultamos el punto final de la API pública de FIVEx (https://fivex.sph.umich.edu/) y recuperamos los resultados en formato JSON. Los resultados se filtraron para retener solo los registros del colon sigmoide (colonsigmoid) y el colon transverso (colontransverse).

Contexto de tumor primario: Los conjuntos de datos de eQTL de TCGA-COAD y adenocarcinoma rectal (READ) se obtuvieron de PancanQTL (^[56]) (https://gonglab.hzau.edu.cn/PancanQTL). Se consideró que un SNP tenía evidencia de eQTL de respaldo si estaba significativamente asociado con la expresión génica en el conjunto de datos de tejido normal del colon o en los conjuntos de datos de tumor primario (P_ < 0,05, FDR < 0,05). Este enfoque inclusivo nos permitió capturar los SNPs con un amplio potencial regulador en diferentes estados fisiológicos.

Construcción y evaluación de 14 modelos de clasificación de sondas de metilación del cáncer

Obtuvimos los datos de metilación TCGA 450k de UCSC Xena (https://xenabrowser.net/datapages/). Después de la clasificación y el análisis estadístico, identificamos 14 tipos de cáncer con al menos 10 muestras normales y tumorales cada uno. Extraímos los datos del nivel de metilación de cg08640619 y cg25982657, junto con las etiquetas de las muestras, de los datos de metilación de estos tipos de cáncer. Utilizando la biblioteca Python scikit (versión 3.12.10), construimos conjuntos de datos de entrenamiento y prueba con la función traintestsplit, estableciendo el parámetro testsize en 0,4 para garantizar una representación equilibrada de las muestras normales y tumorales. Implementamos dos algoritmos de clasificación: LogR y SVM. Para evaluar y comparar su rendimiento durante la fase de selección del modelo, utilizamos la validación cruzada de 10 pliegues en el conjunto de entrenamiento y evaluamos la precisión de la clasificación utilizando la función crossval_score de scikit-learn. Luego, el AUROC, calculado mediante el módulo de métricas de scikit-learn, se utilizó como la métrica principal para evaluar el rendimiento de generalización de los clasificadores seleccionados en datos no vistos. Por último, resumimos y representamos los valores de AUROC de los clasificadores LogR y SVM construidos con cg08640619 y cg25982657 para los 14 tipos de cáncer en gráficos de barras.

Análisis estadístico experimental

El análisis estadístico se realizó utilizando GraphPad Prism (versión 8, GraphPad Software Inc.). Los experimentos se repitieron al menos tres veces, y los resultados se presentan como medias ± DE. Se utilizaron pruebas t o análisis de varianza (ANOVA) para las comparaciones estadísticas, con valores de P < 0,05 considerados estadísticamente significativos.

Para evaluar exhaustivamente el potencial regulador de los SNPs identificados, integramos la evidencia de eQTL tanto de contextos normales como tumorales. Contexto de tejido normal: Se obtuvieron estadísticas resumidas de eQTL en tejido normal del colon a través de la plataforma FIVEx (^[55]) (https://fivex.sph.umich.edu/), que agrega y proporciona datos procesados del Catálogo de eQTL del Instituto Europeo de Bioinformática (EBI) (www.ebi.ac.uk/eqtl/). Para cada variante objetivo, consultamos el punto final de la API pública de FIVEx (https://fivex.sph.umich.edu/) y recuperamos los resultados en formato JSON. Los resultados se filtraron para retener solo los registros del colon sigmoide (colonsigmoid) y el colon transverso (colontransverse).

Contexto de tumor primario: Los conjuntos de datos de eQTL de TCGA-COAD y adenocarcinoma rectal (READ) se obtuvieron de PancanQTL (^[56]) (https://gonglab.hzau.edu.cn/PancanQTL). Se consideró que un SNP tenía evidencia de eQTL de respaldo si estaba significativamente asociado con la expresión génica en el conjunto de datos de tejido normal del colon o en los conjuntos de datos de tumor primario (P_ < 0,05, FDR < 0,05). Este enfoque inclusivo nos permitió capturar SNPs con un amplio potencial regulador en diferentes estados fisiológicos.

Construcción y evaluación de 14 modelos de clasificación de sondas de metilación de cáncer

Obtuvimos los datos de metilación TCGA 450k de UCSC Xena (https://xenabrowser.net/datapages/). Después de la clasificación y el análisis estadístico, identificamos 14 tipos de cáncer con al menos 10 muestras normales y tumorales cada uno. Extraímos los datos del nivel de metilación de cg08640619 y cg25982657, junto con las etiquetas de las muestras, de los datos de metilación de estos tipos de cáncer. Utilizando la biblioteca Python scikit (versión 3.12.10), construimos conjuntos de datos de entrenamiento y prueba con la función traintestsplit, estableciendo el parámetro testsize en 0,4 para garantizar una representación equilibrada de las muestras normales y tumorales. Implementamos dos algoritmos de clasificación: LogR y SVM. Para evaluar y comparar su rendimiento durante la fase de selección del modelo, utilizamos la validación cruzada de 10 pliegues en el conjunto de entrenamiento y evaluamos la precisión de la clasificación utilizando la función crossval_score de scikit-learn. Luego, el AUROC, calculado a través del módulo de métricas de scikit-learn, se utilizó como la métrica principal para evaluar el rendimiento de generalización de los clasificadores seleccionados en datos no vistos. Por último, resumimos y representamos gráficamente los valores de AUROC de los clasificadores LogR y SVM construidos con cg08640619 y cg25982657 para los 14 tipos de cáncer en gráficos de barras.

Análisis estadístico experimental

El análisis estadístico se realizó utilizando GraphPad Prism (versión 8, GraphPad Software Inc.). Los experimentos se repitieron al menos tres veces, y los resultados se presentan como medias ± DE. Se utilizaron pruebas t o análisis de varianza (ANOVA) para las comparaciones estadísticas, y se consideraron estadísticamente significativos los valores de P < 0,05.