El análisis transcriptómico de células individuales revela la heterogeneidad y la dinámica del agotamiento de las células T citotóxicas CD27(+) en el microambiente tumoral del cáncer colorrectal.

El cáncer colorrectal (CCR) representa una de las neoplasias más prevalentes a nivel mundial, ocupando el tercer lugar entre los cánceres más comúnmente diagnosticados y siendo la segunda causa principal de muerte relacionada con el cáncer a nivel mundial ([Sung et al., 2021]). A pesar de los avances en la resección quirúrgica, la quimioterapia y las terapias dirigidas, la tasa de supervivencia a 5 años para el CCR metastásico sigue siendo inferior al 15%, lo que destaca la urgente necesidad de nuevos enfoques terapéuticos ([Sung et al., 2021]). Los recientes avances en la inmunoterapia contra el cáncer, en particular el bloqueo de puntos de control inmunitarios dirigido a las vías PD-1/PD-L1 y CTLA-4, han revolucionado los paradigmas del tratamiento del cáncer y han demostrado una notable eficacia clínica en varios tumores sólidos ([Ribas y Wolchok, 2018]). Sin embargo, las tasas de respuesta a la inmunoterapia en pacientes con CCR siguen siendo muy variables, y solo los tumores con alta inestabilidad de microsatélites (MSI-H) muestran un beneficio significativo, mientras que los tumores con estabilidad de microsatélites (MSS), que constituyen aproximadamente el 85% de los casos de CCR, siguen siendo en gran medida resistentes a las estrategias inmunoterapéuticas actuales ([Le et al., 2017]).

El microambiente tumoral (MET) desempeña un papel fundamental en la determinación de las respuestas terapéuticas y los resultados clínicos en el CCR ([Kather y Halama, 2019]). Los linfocitos T citotóxicos, en particular las células T CD8+, representan las principales células efectoras responsables de la inmunidad antitumoral a través del reconocimiento directo y la eliminación de las células malignas ([Zhang et al., 2018]). Entre las poblaciones de células T CD8+, las células T CD27+ han surgido como un subconjunto crítico con distintas propiedades funcionales y significado pronóstico en varios cánceres ([Hendriks et al., 2000]). CD27, un miembro de la superfamilia de receptores del factor de necrosis tumoral, actúa como una importante molécula costimuladora que regula la activación, la proliferación y la supervivencia de las células T, y la señalización de CD27 promueve respuestas inmunitarias antitumorales sostenidas ([van de Ven y Borst, 2015]). Estudios previos han demostrado que la expresión de CD27 se correlaciona con una mayor funcionalidad de las células T y resultados clínicos favorables en múltiples tipos de cáncer, lo que sugiere su potencial como biomarcador pronóstico y diana terapéutica ([Claus et al., 2012]).

A pesar de la importancia reconocida de las células T CD27+ en la inmunidad antitumoral, un desafío fundamental en la inmunoterapia contra el cáncer radica en el fenómeno del agotamiento de las células T, un estado disfuncional caracterizado por la pérdida progresiva de las funciones efectoras, la expresión sostenida de múltiples receptores inhibitorios y la capacidad proliferativa alterada ([Wherry y Kurachi, 2015]). El agotamiento de las células T se desarrolla a través de la estimulación antigénica crónica dentro del MET inmunosupresor y representa una barrera importante para una inmunidad antitumoral eficaz ([McLane et al., 2019]). Las células T agotadas exhiben un programa transcripcional distinto orquestado por factores de transcripción clave, incluidos TOX, TCF1 y EOMES, que impulsan la expresión de moléculas de puntos de control inmunitarios como PD-1, LAG-3, TIM-3 y CTLA-4 ([Khan et al., 2019]). Comprender los mecanismos moleculares y la heterogeneidad celular que subyacen al agotamiento de las células T es esencial para desarrollar estrategias para revitalizar las respuestas inmunitarias antitumorales y mejorar la eficacia de la inmunoterapia.

Las tecnologías de secuenciación de ARN de una sola célula (scRNA-seq) han revolucionado nuestra comprensión de la heterogeneidad celular dentro del MET, lo que permite una resolución sin precedentes de distintos estados celulares, trayectorias de diferenciación y programas funcionales a nivel de célula individual ([Papalexi y Satija, 2018]). Estudios recientes de scRNA-seq en varios cánceres han revelado una notable heterogeneidad dentro de las poblaciones de células T agotadas, identificando distintos subconjuntos agotados con perfiles transcripcionales, capacidades funcionales y capacidad de respuesta terapéutica diferenciales ([Thommen et al., 2018]). Además, los algoritmos de inferencia de trayectorias aplicados a los datos de scRNA-seq han dilucidado los procesos dinámicos de diferenciación de las células T y la progresión del agotamiento, revelando puntos de decisión críticos y reguladores moleculares que gobiernan la determinación del destino celular ([Qiu et al., 2017]). Sin embargo, la caracterización integral a nivel de una sola célula de la heterogeneidad de las células T citotóxicas CD27+, la dinámica del agotamiento y la organización espacial específicamente dentro del MET del CCR sigue siendo limitada.

En este estudio, realizamos un análisis transcriptómico integral de una sola célula de las células T citotóxicas CD27+ aisladas de los tejidos tumorales del CCR y de los tejidos adyacentes normales correspondientes para analizar la heterogeneidad celular, identificar distintos estados funcionales, reconstruir las trayectorias de diferenciación y dilucidar los mecanismos moleculares que subyacen al agotamiento de las células T. A través de enfoques computacionales integradores que incluyen la reducción de dimensionalidad, la inferencia de trayectorias, el análisis de redes de regulación génica y el mapeo de interacciones célula-célula, nuestro objetivo fue proporcionar un atlas de alta resolución de los estados de las células T CD27+ en el CCR e identificar posibles dianas terapéuticas para revertir el agotamiento de las células T y mejorar la inmunidad antitumoral.

Métodos

Adquisición y preprocesamiento de datos

Los datos de secuenciación de ARN de una sola célula se obtuvieron de la base de datos Gene Expression Omnibus (GEO) con el número de acceso GSE144735. El conjunto de datos comprendía poblaciones de células T citotóxicas CD27+ aisladas de tejidos tumorales de cáncer colorrectal (CCR) y tejidos adyacentes normales correspondientes. Los datos de secuenciación brutos se sometieron a un filtrado de control de calidad para eliminar las células y los dobletes de baja calidad. Las células con menos de 200 genes detectados, más del 10% de contenido de genes mitocondriales o recuentos extremos de identificadores moleculares únicos (UMI) se excluyeron del análisis posterior. Las matrices de expresión génica se normalizaron utilizando procedimientos estándar de normalización logarítmica, y se identificaron los genes altamente variables para los análisis posteriores de reducción de dimensionalidad.

Reducción de dimensionalidad y análisis de agrupamiento

Se emplearon múltiples técnicas complementarias de reducción de dimensionalidad para capturar la heterogeneidad multidimensional de las poblaciones de células T CD27+. Se aplicó la aproximación y proyección de múltiples uniformes (UMAP) para visualizar los grupos celulares discretos, mientras que la incrustación de vecinos estocásticos distribuidos en t (t-SNE) se utilizó para capturar los gradientes continuos de los estados celulares, en particular el gradiente de agotamiento desde fenotipos funcionales hasta fenotipos terminalmente agotados. Se realizó un análisis de componentes principales (PCA) para identificar los principales impulsores moleculares de la heterogeneidad de las células T, y el análisis de contribución de la varianza determinó el número óptimo de componentes principales. Se seleccionaron los tres primeros componentes principales, que capturan aproximadamente el 10% de la varianza acumulada, basándose en el análisis del codo de la gráfica de dispersión en PC8. Se generaron proyecciones de PCA tridimensionales para validar la naturaleza multidimensional de la heterogeneidad celular. Se aplicaron algoritmos de agrupamiento basados en gráficos para identificar poblaciones celulares discretas, y los límites de los grupos se delimitaron utilizando el análisis de envolvente convexa. Además, se emplearon diseños de gráficos dirigidos por fuerzas y el mapeo de difusión para ilustrar las relaciones interconectadas entre los tipos de células y para capturar el paisaje de la función citotóxica a lo largo de gradientes continuos.

Inferencia de trayectorias y análisis de pseudotiempo

Para reconstruir las trayectorias de diferenciación de las células T y la progresión del agotamiento, se implementaron múltiples algoritmos de inferencia de trayectorias. Se realizó un análisis de pseudotiempo para ordenar las células a lo largo del continuo de agotamiento, desde estados terminalmente agotados hasta estados efectores. La construcción del árbol de trayectorias basada en Monocle identificó los eventos de ramificación que conducen a distintos destinos de diferenciación, incluidos los caminos efectores, de memoria y agotados. Se ajustaron curvas principales de Slingshot a través del espacio de componentes principales para capturar trayectorias lineales suaves para cada ruta de diferenciación. El análisis de conectividad basado en la abstracción de gráficos de partición (PAGA) construyó redes integrales de transición de estados, revelando relaciones jerárquicas entre los estados ingenuo, activado, de memoria, efector y agotado. Las probabilidades de transición entre los estados funcionales se cuantificaron utilizando el análisis del flujo de transición de estados celulares. Se realizó un análisis de velocidad de ARN para inferir la dinámica transcripcional direccional y predecir los estados celulares futuros, y las líneas de flujo de velocidad indicaron las trayectorias de diferenciación predominantes y las decisiones sobre el destino celular.

Expresión génica diferencial y caracterización funcional

Se realizó un análisis de expresión génica diferencial para identificar los genes marcadores que distinguen los subconjuntos de células T CD27+ y para comparar la expresión génica entre los tejidos tumorales y los tejidos adyacentes normales. La significación estadística se evaluó utilizando pruebas de suma de rangos de Wilcoxon con corrección de Bonferroni para pruebas múltiples. Se consideraron significativamente diferentes los genes con un cambio de logaritmo de 2 mayor que 0,5 y un valor de p ajustado menor que 0,05. Los resultados se visualizaron utilizando gráficos de volcanes, gráficos de MA y gráficos de cascada para resaltar los principales marcadores de agotamiento (PDCD1, CTLA4, HAVCR2, LAG3) y las moléculas efectoras. Los patrones de expresión génica en los grupos se mostraron utilizando mapas de calor con agrupamiento jerárquico para revelar las firmas moleculares de las poblaciones terminales agotadas, efectoras y similares a la memoria. Las matrices de gráficos de puntos perfilaron sistemáticamente la expresión de marcadores específicos del tipo de célula, con el tamaño del punto que indica el porcentaje de células que expresan y el color que representa los niveles de expresión promedio. Se generaron gráficos de violín, gráficos de caja y gráficos de cresta para comparar las distribuciones de expresión de los genes críticos (CD27, CD8A, PRF1, CD3E, GZMB, TOX, TCF7) en los estados funcionales y entre los tipos de tejido.

Puntuación funcional y análisis de enriquecimiento de vías

Se cuantificaron las características funcionales de los subconjuntos de células T utilizando enfoques de puntuación basados en conjuntos de genes. Se calcularon las puntuaciones de agotamiento mediante la agregación de la expresión de los genes de la firma de agotamiento establecida, mientras que las puntuaciones de citotoxicidad integraron la expresión de las moléculas efectoras citotóxicas, incluidas las granzimas y la perforina. Se calcularon de manera similar dimensiones funcionales adicionales, que incluyen las firmas de activación, proliferación y memoria, para un perfil completo en gráficos de radar en los subtipos de células T CD8+. Se realizó un análisis de enriquecimiento de conjuntos de genes (GSEA) para identificar las vías y los procesos biológicos enriquecidos en los distintos estados de las células T, y se determinó la puntuación de enriquecimiento normalizada (NES) y la significación estadística mediante pruebas de permutación. El análisis de enriquecimiento de vías en múltiples procesos biológicos, que incluyen la señalización de citocinas, la reprogramación metabólica y la regulación del ciclo celular, se visualizó utilizando gráficos de burbujas con puntuaciones de enriquecimiento y valores de significación.

Análisis de la red de regulación génica y la interacción célula-célula

Se construyeron redes reguladoras de genes para mapear las interacciones entre factores de transcripción y genes diana que gobiernan el agotamiento de las células T y sus funciones efectoras. La inferencia de la red reguladora de genes (GRN) se realizó utilizando SCENIC (Inferencia y agrupación de redes reguladoras a nivel de célula única), implementado a través del pipeline pySCENIC (pySCENIC v0.12.1). Se llevó a cabo un análisis de los motivos de unión de los factores de transcripción utilizando la base de datos de motivos cisTarget/JASPAR (JASPAR 2022) para identificar motivos reguladores enriquecidos e inferir las relaciones entre los factores de transcripción y los genes diana. Se aplicó un umbral de correlación de Pearson de r ≥ 0,3 para la construcción del módulo de coexpresión, y se realizó la poda del regulón utilizando el paso RcisTarget con un umbral de puntuación de enriquecimiento normalizado (NES) de 3,0. Las aristas activadoras se definieron como correlaciones positivas entre la expresión del factor de transcripción y el gen diana, mientras que las aristas represoras se definieron como correlaciones negativas con evidencia de motivos de respaldo. Los gráficos de red representaron los nodos reguladores clave que controlan las transiciones de efector a agotado, con aristas que representan interacciones activadoras o represoras. Se realizó un análisis de la interacción ligando-receptor para mapear las vías de comunicación entre células dentro del microambiente tumoral del CCRC. Se identificaron posibles interacciones entre las células T, las células mieloides y otras poblaciones inmunitarias basándose en la coexpresión de pares de ligandos y receptores cognados, incluidas las interacciones de los puntos de control inmunitario (PD-1/PD-L1, CTLA-4/CD80/CD86) y las vías de señalización de las citocinas. Se cuantificaron y visualizaron las intensidades de interacción utilizando diagramas de red y matrices de comunicación.

Análisis espacial y composicional

Se cuantificaron las proporciones de tipos de células en todos los grupos de muestras y se compararon entre los tejidos tumorales y los tejidos adyacentes normales utilizando gráficos de barras apiladas y análisis composicional. La agrupación jerárquica de los patrones de firmas celulares en los diferentes grupos de muestras reveló una heterogeneidad específica del paciente en la composición de las células T a través de la visualización de clustermap. Se modelaron los patrones de distribución espacial utilizando un enfoque de simulación computacional basado en los hallazgos de transcriptómica espacial publicados en CCRC, ya que no se dispuso de datos de transcriptómica espacial o inmunofluorescencia multiplex directa en el conjunto de datos GSE144735. Este modelo espacial simulado se presenta explícitamente como una inferencia computacional consistente con los patrones espaciales establecidos del CCRC y no representa una medición espacial empírica directa. Se realizó un análisis del clonotipo del receptor de células T (TCR) para evaluar los patrones de expansión clonal y la diversidad, identificando los principales clonotipos expandidos y visualizando su distribución en los diferentes grupos de muestras utilizando gráficos de barras y gráficos circulares.

Dinámica de la expresión génica temporal

Se analizaron las dinámicas de la expresión génica a lo largo de las trayectorias de agotamiento ordenando las células según el pseudotiempo y examinando los patrones de regulación temporal. Los mapas de calor de la cascada de expresión génica revelaron distintas fases temporales, con genes efectores tempranos (CD69, PRF1) que transicionaron a la expresión intermedia de puntos de control (PDCD1, CTLA4) y marcadores de agotamiento terminal (HAVCR2, TOX, LAG3). El análisis de correlación cuantificó las relaciones entre las puntuaciones de agotamiento y la citotoxicidad a lo largo de la progresión del pseudotiempo. Los gráficos de doble eje demostraron correlaciones inversas entre las medidas funcionales y los estados de agotamiento.

Experimentos de cocultivo de líneas celulares

Para validar los hallazgos computacionales del análisis de secuenciación de ARN de célula única, se realizaron experimentos de cocultivo in vitro utilizando cuatro líneas celulares de cáncer colorrectal con diferentes características moleculares: HCT116 y RKO (inestabilidad de microsatélites alta, MSI-H) y HCT15 y SW480 (inestabilidad de microsatélites estable, MSS). Las células cancerosas se cultivaron en medio RPMI-1640 suplementado con un 10% de suero fetal bovino y un 1% de penicilina-estreptomicina a 37 °C en una atmósfera humidificada que contenía un 5% de CO2. Se obtuvieron células mononucleares de sangre periférica humana (PBMC) comercialmente de la American Type Culture Collection (ATCC; número de catálogo PCS-800-011) y se utilizaron de acuerdo con las directrices del proveedor. Las PBMC se descongelaron y se recuperaron en medio RPMI-1640 completo durante 2 horas antes de la activación. Las PBMC se activaron con anticuerpos anti-CD3/CD28 (1 μg/mL cada uno) durante 24 horas antes del cocultivo. Las células cancerosas se sembraron a una densidad de 1 × 105 células por pocillo en placas de 6 pocillos y se dejaron adherir durante la noche. A continuación, se añadieron las PBMC activadas a una relación efector-diana de 10:1, y los cocultivos se mantuvieron durante 48 horas en condiciones de cultivo estándar.

Tras el período de cocultivo, las PBMC se recogieron cuidadosamente y se extrajo el ARN total utilizando el reactivo TRIzol según el protocolo del fabricante. Se evaluó la concentración y la pureza del ARN mediante espectrofotometría NanoDrop, y se consideraron aceptables las relaciones A260/A280 entre 1,8 y 2,0. Se sintetizó ADN complementario (ADNc) a partir de 1 μg de ARN total utilizando un kit de transcripción inversa con cebadores oligo (dT). Se realizó una PCR cuantitativa de transcripción inversa (RT-qPCR) en un sistema de PCR en tiempo real utilizando la mezcla maestra SYBR Green para cuantificar la expresión de los marcadores de agotamiento de las células T, incluidos PDCD1 (que codifica PD-1), LAG3, HAVCR2 (que codifica TIM-3) y CD27. Se utilizaron las siguientes secuencias de cebadores para la amplificación de RT-qPCR: PDCD1 directo: 5′-CCA​GGA​TGG​TTC​TTA​GAC​TCC​C-3′, reverso: 5′-TTT​AGC​ACG​AAG​CTC​TCC​GAT-3'; LAG3 directo: 5′-CTG​GAT​GGT​TTC​TGC​GAG​AC-3′, reverso: 5′-CCA​GTC​ATC​TTC​TGG​TGG​GA-3'; HAVCR2 directo: 5′-TTC​CCG​ACA​AAG​GCT​TGA​CT-3′, reverso: 5′-GCT​GGG​TTC​TTG​GTC​AGC​AT-3'; CD27 directo: 5′-CTG​CAA​CCA​GAG​AAA​GGG​TAG​A-3′, reverso: 5′-GCT​TTG​CAT​CTC​TTG​GTT​GAG​A-3'. Se utilizó GAPDH como gen de referencia interno con los cebadores: directo: 5′-GTC​TCC​TCT​GAC​TTC​AAC​AGC​G-3′, reverso: 5′-ACC​ACC​CTG​TTG​CTG​TAG​CCA​A-3'. Todos los cebadores se diseñaron para abarcar uniones exón-exón para evitar la amplificación del ADN genómico y se validaron para la especificidad mediante el análisis de la curva de fusión y la electroforesis en gel de agarosa.

Las reacciones de RT-qPCR se realizaron en triplicado técnico con las siguientes condiciones de ciclado: desnaturalización inicial a 95 °C durante 10 minutos, seguida de 40 ciclos de desnaturalización a 95 °C durante 15 segundos y hibridación/extensión a 60 °C durante 60 segundos. Se calculó la expresión génica relativa utilizando el método 2-ΔΔCt, y se determinaron los cambios en relación con los valores de referencia del monocultivo, donde las células cancerosas o las PBMC se cultivaron de forma independiente sin compañeros de cocultivo. El valor de ΔCt se calculó restando el valor de Ct de GAPDH del valor de Ct del gen diana para cada muestra, y ΔΔCt se calculó restando el ΔCt del control de monocultivo del ΔCt de la muestra de cocultivo. Los datos representan los valores medios de tres réplicas biológicas independientes, cada una con triplicados técnicos, y los resultados se expresaron como cambios en relación con el valor de referencia del monocultivo (establecido en 1,0). Se realizaron comparaciones estadísticas entre las condiciones de cocultivo MSI-H (HCT116 y RKO) y MSS (HCT15 y SW480) utilizando pruebas t de Student de dos colas para comparaciones por pares y ANOVA de una vía con la prueba de Tukey post hoc para comparaciones de múltiples grupos, y se consideraron estadísticamente significativos los valores de p inferiores a 0,05.

Correlación clínica y análisis de supervivencia

Se realizó un análisis de supervivencia para evaluar de forma exploratoria la posible relevancia pronóstica de las firmas de agotamiento de las células T en un conjunto de datos disponible, y los resultados tenían como objetivo ser únicamente generadores de hipótesis. Los pacientes se estratificaron en grupos de agotamiento alto y bajo en función de las puntuaciones de agotamiento medianas, y se generaron curvas de supervivencia de Kaplan-Meier. Se calcularon las razones de riesgo y los valores de p utilizando modelos de riesgos proporcionales de Cox. Se realizó un cálculo de potencia post hoc para evaluar el tamaño de muestra necesario para detectar el tamaño del efecto observado (HR = 0,505) con una potencia del 80% a un alfa bilateral de 0,05, contextualizando la probable falta de potencia del grupo disponible para este punto final. Se comparó la expresión del objetivo terapéutico entre las poblaciones de células T agotadas y efectoras para identificar las moléculas de puntos de control que mostraban una sobreexpresión significativa en las células agotadas, lo que proporcionaba información para la posible priorización de los objetivos de la inmunoterapia.

Análisis estadístico y visualización

Todos los análisis estadísticos se realizaron utilizando R (versión 4.0 o superior) con Seurat, Scanpy y paquetes de bioinformática relacionados. Se utilizaron pruebas estadísticas de dos colas en todo momento, y se consideraron estadísticamente significativos los valores de p inferiores a 0,05, a menos que se especifique lo contrario. Se aplicó la corrección de pruebas múltiples utilizando el método de Bonferroni o Benjamini-Hochberg según corresponda. La visualización de datos se generó utilizando ggplot2, ComplexHeatmap y funciones de trazado personalizadas, y las figuras se diseñaron para integrar múltiples perspectivas analíticas, incluido el análisis composicional, la inferencia de trayectoria, la elaboración de perfiles funcionales y las correlaciones clínicas.

Adquisición y preprocesamiento de datos

Los datos de secuenciación de ARN de célula única se obtuvieron de la base de datos Gene Expression Omnibus (GEO) con el número de acceso GSE144735. El conjunto de datos comprendía poblaciones de células T citotóxicas CD27+ aisladas de tejidos tumorales de cáncer colorrectal (CCRC) y tejidos adyacentes normales. Los datos de secuenciación brutos se sometieron a un filtrado de control de calidad para eliminar las células y los dobletes de baja calidad. Se excluyeron del análisis posterior las células con menos de 200 genes detectados, más del 10% de contenido de genes mitocondriales o conteos extremos de identificadores moleculares únicos (UMI). Las matrices de expresión génica se normalizaron utilizando procedimientos de normalización logarítmica estándar, y se identificaron los genes altamente variables para los análisis posteriores de reducción de dimensionalidad.

Reducción de dimensionalidad y análisis de agrupación

Se emplearon múltiples técnicas complementarias de reducción de dimensionalidad para capturar la heterogeneidad multidimensional de las poblaciones de células T CD27+. Se aplicó la aproximación y proyección de múltiples uniformes (UMAP) para visualizar los grupos celulares discretos, mientras que la incrustación de vecinos estocásticos distribuidos en t (t-SNE) se utilizó para capturar los gradientes continuos de los estados celulares, en particular el gradiente de agotamiento desde fenotipos funcionales hasta fenotipos terminalmente agotados. Se realizó un análisis de componentes principales (PCA) para identificar los principales impulsores moleculares de la heterogeneidad de las células T, y el análisis de la contribución de la varianza determinó el número óptimo de componentes principales. Se seleccionaron los tres primeros componentes principales, que capturaron aproximadamente el 10% de la varianza acumulada, basándose en el análisis del codo de la gráfica de dispersión en PC8. Se generaron proyecciones de PCA tridimensionales para validar la naturaleza multidimensional de la heterogeneidad celular. Se aplicaron algoritmos de agrupación basados en gráficos para identificar poblaciones celulares discretas, y los límites de los grupos se delimitaron utilizando el análisis de la envolvente convexa. Además, se utilizaron diseños de gráficos dirigidos por fuerzas y el mapeo de difusión para ilustrar las relaciones interconectadas entre los tipos de células y para capturar el panorama de la función citotóxica a lo largo de gradientes continuos.

Inferencia de trayectoria y análisis de pseudotiempo

Para reconstruir las trayectorias de diferenciación de las células T y la progresión del agotamiento, se implementaron múltiples algoritmos de inferencia de trayectorias. Se realizó un análisis de pseudotiempo para ordenar las células a lo largo del continuo de agotamiento, desde estados de agotamiento terminal hasta estados efectores. La construcción de un árbol de trayectorias basado en Monocle identificó eventos de ramificación que conducen a destinos de diferenciación distintos, incluidos los linajes efector, de memoria y de agotamiento. Se ajustaron curvas principales de Slingshot a través del espacio de componentes principales para capturar trayectorias lineales suaves para cada ruta de diferenciación. El análisis de conectividad basado en la abstracción de grafos por partición (PAGA) construyó redes integrales de transición de estados, revelando relaciones jerárquicas entre los estados naive, activado, de memoria, efector y de agotamiento. Las probabilidades de transición entre los estados funcionales se cuantificaron utilizando el análisis del flujo de transición de estados celulares. Se realizó un análisis de la velocidad del ARN para inferir la dinámica transcripcional direccional y predecir los estados celulares futuros, con líneas de flujo de velocidad que indican las trayectorias de diferenciación predominantes y las decisiones sobre el destino celular.

Expresión génica diferencial y caracterización funcional

Se realizó un análisis de expresión génica diferencial para identificar los genes marcadores que distinguen las subpoblaciones de células T CD27+ y para comparar la expresión génica entre los tejidos tumorales y los tejidos adyacentes normales. La significación estadística se evaluó utilizando pruebas de rango-suma de Wilcoxon con corrección de Bonferroni para pruebas múltiples. Los genes con un cambio de pliegue log2 mayor que 0,5 y un valor p ajustado menor que 0,05 se consideraron significativamente expresados de forma diferencial. Los resultados se visualizaron utilizando gráficos de volcán, gráficos MA y gráficos de cascada para resaltar los principales marcadores de agotamiento (PDCD1, CTLA4, HAVCR2, LAG3) y las moléculas efectoras. Los patrones de expresión génica en los grupos se mostraron utilizando mapas de calor con agrupamiento jerárquico para revelar las firmas moleculares de las poblaciones de agotamiento terminal, efectoras y similares a la memoria. Las matrices de gráficos de puntos perfilaron sistemáticamente la expresión de marcadores específicos de cada tipo de célula, con el tamaño del punto que indica el porcentaje de células que expresan y el color que representa los niveles de expresión promedio. Se generaron gráficos de violín, gráficos de caja y gráficos de cresta para comparar las distribuciones de expresión de los genes críticos (CD27, CD8A, PRF1, CD3E, GZMB, TOX, TCF7) en los diferentes estados funcionales y entre los tipos de tejido.

Puntuación funcional y análisis de enriquecimiento de vías

Se cuantificaron las características funcionales de los subconjuntos de células T utilizando enfoques de puntuación basados en conjuntos de genes. Los puntajes de agotamiento se calcularon agregando la expresión de los genes de la firma de agotamiento establecida, mientras que los puntajes de citotoxicidad integraron la expresión de las moléculas efectoras citotóxicas, incluidas las granzimas y la perforina. Se calcularon de manera similar dimensiones funcionales adicionales, que incluyen las firmas de activación, proliferación y memoria, para un perfil completo en gráficos de radar en los subtipos de células T CD8+. Se realizó un análisis de enriquecimiento de conjuntos de genes (GSEA) para identificar las vías y los procesos biológicos enriquecidos en los diferentes estados de las células T, con una puntuación de enriquecimiento normalizada (NES) y una significación estadística determinadas mediante pruebas de permutación. El análisis de enriquecimiento de vías en múltiples procesos biológicos, que incluyen la señalización de citocinas, la reprogramación metabólica y la regulación del ciclo celular, se visualizó utilizando gráficos de burbujas con puntajes de enriquecimiento y valores de significación.

Análisis de la red reguladora génica y la interacción célula-célula

Se construyeron redes reguladoras génicas para mapear las interacciones entre los genes diana y los factores de transcripción que gobiernan el agotamiento de las células T y las funciones efectoras. La inferencia de la GRN se realizó utilizando SCENIC (Inferencia y agrupamiento de la red reguladora de una sola célula), implementada a través de la canalización pySCENIC (pySCENIC v0.12.1). Se realizó un análisis del motivo de unión del factor de transcripción utilizando la base de datos de motivos cisTarget/JASPAR (JASPAR 2022) para identificar los motivos reguladores enriquecidos e inferir las relaciones entre los genes diana y los factores de transcripción. Se aplicó un umbral de correlación de Pearson de r ≥ 0,3 para la construcción del módulo de coexpresión, y se realizó la poda del regulón utilizando el paso RcisTarget con un umbral de puntuación de enriquecimiento normalizada (NES) de 3,0. Los bordes activadores se definieron como correlaciones positivas entre la expresión del factor de transcripción y el gen diana, mientras que los bordes represores se definieron como correlaciones negativas con evidencia de motivos de respaldo. Los gráficos de red representaron los nodos reguladores clave que controlan las transiciones de efector a agotamiento, con bordes que representan las interacciones activadoras o represoras. Se realizó un análisis de la interacción ligando-receptor para mapear las vías de comunicación célula-célula dentro del microambiente tumoral del CRC. Se identificaron interacciones potenciales entre las células T, las células mieloides y otras poblaciones inmunitarias en función de la coexpresión de pares de ligandos y receptores cognados, incluidas las interacciones de los puntos de control inmunitario (PD-1/PD-L1, CTLA-4/CD80/CD86) y las vías de señalización de citocinas. Se cuantificaron y visualizaron las intensidades de interacción utilizando diagramas de red y matrices de comunicación.

Análisis espacial y composicional

Se cuantificaron las proporciones de tipos de células en las muestras y se compararon entre los tejidos tumorales y los tejidos adyacentes normales utilizando gráficos de barras apiladas y análisis composicional. El agrupamiento jerárquico de los patrones de firma celular en las muestras individuales reveló la heterogeneidad específica del paciente en la composición de las células T a través de la visualización de clustermap. Se modelaron los patrones de distribución espacial utilizando un enfoque de simulación computacional basado en los hallazgos de transcriptómica espacial publicados en CRC, ya que no estaban disponibles datos directos de transcriptómica espacial o inmunofluorescencia multiplex en el conjunto de datos GSE144735. Este modelo espacial simulado se presenta explícitamente como una inferencia computacional consistente con los patrones espaciales establecidos del CRC y no representa una medición espacial empírica directa. Se realizó un análisis del clonotipo del receptor de células T (TCR) para evaluar los patrones de expansión clonal y la diversidad, con los principales clonotipos expandidos identificados y su distribución en las muestras visualizada utilizando gráficos de barras y gráficos circulares.

Dinámica de la expresión génica temporal

Se analizaron las dinámicas de la expresión génica a lo largo de las trayectorias de agotamiento ordenando las células según el pseudotiempo y examinando los patrones de regulación temporal. Los mapas de calor de la cascada de expresión génica revelaron distintas fases temporales, con los primeros genes efectores (CD69, PRF1) que hacen la transición a la expresión intermedia de puntos de control (PDCD1, CTLA4) y los marcadores de agotamiento terminal (HAVCR2, TOX, LAG3). El análisis de correlación cuantificó las relaciones entre los puntajes de agotamiento y la citotoxicidad a lo largo de la progresión del pseudotiempo. Los gráficos de doble eje demostraron correlaciones inversas entre las medidas funcionales y los estados de agotamiento.

Experimentos de cocultivo de líneas celulares

Para validar los hallazgos computacionales del análisis de secuenciación de ARN de una sola célula, se realizaron experimentos de cocultivo in vitro utilizando cuatro líneas celulares de cáncer colorrectal con diferentes características moleculares: HCT116 y RKO (inestabilidad de microsatélites alta, MSI-H) y HCT15 y SW480 (inestabilidad de microsatélites estable, MSS). Las células cancerosas se cultivaron en medio RPMI-1640 suplementado con un 10% de suero fetal bovino y un 1% de penicilina-estreptomicina a 37 °C en una atmósfera humidificada que contiene un 5% de CO2. Las células mononucleares de sangre periférica humana (PBMC) se obtuvieron comercialmente de la American Type Culture Collection (ATCC; número de catálogo PCS-800-011) y se utilizaron de acuerdo con las pautas del proveedor. Las PBMC se descongelaron y se recuperaron en medio completo RPMI-1640 durante 2 horas antes de la activación. Las PBMC se activaron con anticuerpos anti-CD3/CD28 (1 μg/mL cada uno) durante 24 horas antes del cocultivo. Las células cancerosas se sembraron a una densidad de 1 × 105 células por pocillo en placas de 6 pocillos y se dejaron adherir durante la noche. Luego, se agregaron las PBMC activadas a una relación efector-diana de 10:1, y los cocultivos se mantuvieron durante 48 horas en condiciones de cultivo estándar.

Después del período de cocultivo, las PBMC se recolectaron cuidadosamente y se extrajo el ARN total utilizando el reactivo TRIzol según el protocolo del fabricante. La concentración y la pureza del ARN se evaluaron mediante espectrofotometría NanoDrop, y las relaciones A260/A280 entre 1,8 y 2,0 se consideraron aceptables. El ADN complementario (ADNc) se sintetizó a partir de 1 μg de ARN total utilizando un kit de transcripción inversa con cebadores de oligo (dT). Se realizó una reacción en cadena de la polimerasa cuantitativa con transcripción inversa (RT-qPCR) en un sistema de PCR en tiempo real utilizando la mezcla maestra SYBR Green para cuantificar la expresión de los marcadores de agotamiento de las células T, incluidos PDCD1 (que codifica PD-1), LAG3, HAVCR2 (que codifica TIM-3) y CD27. Se utilizaron las siguientes secuencias de cebadores para la amplificación de RT-qPCR: PDCD1 directo: 5′-CCA​GGA​TGG​TTC​TTA​GAC​TCC​C-3′, reverso: 5′-TTT​AGC​ACG​AAG​CTC​TCC​GAT-3'; LAG3 directo: 5′-CTG​GAT​GGT​TTC​TGC​GAG​AC-3′, reverso: 5′-CCA​GTC​ATC​TTC​TGG​TGG​GA-3'; HAVCR2 directo: 5′-TTC​CCG​ACA​AAG​GCT​TGA​CT-3′, reverso: 5′-GCT​GGG​TTC​TTG​GTC​AGC​AT-3'; CD27 directo: 5′-CTG​CAA​CCA​GAG​AAA​GGG​TAG​A-3′, reverso: 5′-GCT​TTG​CAT​CTC​TTG​GTT​GAG​A-3'. GAPDH se utilizó como gen de referencia interno con los cebadores: directo: 5′-GTC​TCC​TCT​GAC​TTC​AAC​AGC​G-3′, reverso: 5′-ACC​ACC​CTG​TTG​CTG​TAG​CCA​A-3'. Todos los cebadores se diseñaron para abarcar las uniones exón-exón para evitar la amplificación del ADN genómico y se validaron para la especificidad mediante el análisis de la curva de fusión y la electroforesis en gel de agarosa.

Las reacciones de RT-qPCR se realizaron en triplicado técnico con las siguientes condiciones de ciclado: desnaturalización inicial a 95 °C durante 10 minutos, seguida de 40 ciclos de desnaturalización a 95 °C durante 15 segundos y hibridación/extensión a 60 °C durante 60 segundos. La expresión génica relativa se calculó utilizando el método 2-ΔΔCt, con los cambios de pliegue determinados en relación con las líneas de base del monocultivo donde las células cancerosas o las PBMC se cultivaron de forma independiente sin socios de cocultivo. El valor de ΔCt se calculó restando el valor de Ct de GAPDH del valor de Ct del gen diana para cada muestra, y ΔΔCt se calculó restando el ΔCt del control de monocultivo del ΔCt de la muestra de cocultivo. Los datos representan los valores medios de tres réplicas biológicas independientes, cada una con triplicados técnicos, y los resultados se expresaron como un cambio de pliegue en relación con la línea de base del monocultivo (establecida en 1,0). Las comparaciones estadísticas entre las condiciones de cocultivo MSI-H (HCT116 y RKO) y MSS (HCT15 y SW480) se realizaron utilizando pruebas t de Student de dos colas para comparaciones por pares y ANOVA de una vía con la prueba post hoc de Tukey para comparaciones de múltiples grupos, con valores p menores que 0,05 considerados estadísticamente significativos.

Correlación clínica y análisis de supervivencia

Se realizó un análisis de supervivencia para evaluar de forma exploratoria la posible relevancia pronóstica de las firmas de agotamiento de las células T en un conjunto de datos disponible, y los resultados tenían la intención de ser solo generadores de hipótesis. Los pacientes se estratificaron en grupos de agotamiento alto y bajo en función de los puntajes de agotamiento medianos, y se generaron curvas de supervivencia de Kaplan-Meier. Los riesgos y los valores p se calcularon utilizando modelos de riesgos proporcionales de Cox. Se realizó un cálculo de potencia post hoc para evaluar el tamaño de muestra necesario para detectar el tamaño del efecto observado (HR = 0,505) con una potencia del 80% a un alfa bilateral de 0,05, contextualizando la probable naturaleza poco potente de la cohorte disponible para este punto final. Se comparó la expresión del objetivo terapéutico entre las poblaciones de células T agotadas y efectoras para identificar las moléculas de puntos de control que muestran una sobreexpresión significativa en las células agotadas, lo que informa la posible priorización de los objetivos de la inmunoterapia.

Análisis estadístico y visualización

Todos los análisis estadísticos se realizaron utilizando R (versión 4.0 o superior) con Seurat, Scanpy y paquetes de bioinformática relacionados. Se emplearon pruebas estadísticas de dos colas en todos los casos, y se consideraron estadísticamente significativos los valores p menores de 0,05, a menos que se especifique lo contrario. Se aplicó una corrección para pruebas múltiples utilizando el método de Bonferroni o Benjamini-Hochberg, según corresponda. La visualización de datos se generó utilizando ggplot2, ComplexHeatmap y funciones de trazado personalizadas, y las figuras se diseñaron para integrar múltiples perspectivas analíticas, incluido el análisis composicional, la inferencia de trayectorias, la elaboración de perfiles funcionales y las correlaciones clínicas.

Resultados

Arquitectura celular y análisis de reducción de dimensionalidad de las poblaciones de células T CD27+

La caracterización inicial revela patrones de agrupamiento distintos de las células T citotóxicas CD27+ dentro del microambiente tumoral del cáncer colorrectal. El análisis UMAP (Figura 1A) identifica múltiples grupos celulares discretos con diversos estados fenotípicos, lo que demuestra una heterogeneidad sustancial dentro del compartimento de células T CD27+. La visualización t-SNE (Figura 1B) captura un gradiente de agotamiento, con células que pasan de estados funcionales a estados de agotamiento terminal, como lo indica el gradiente de color que representa las puntuaciones de agotamiento. La Figura 1C presenta un análisis de contribución génica basado en PCA, que revela los principales impulsores moleculares de la heterogeneidad de las células T. La distribución de la densidad celular (Figura 1D) destaca los patrones de enriquecimiento en todo el microambiente tisular, mientras que el análisis comparativo entre el tejido tumoral y el tejido adyacente normal (Figura 1E) demuestra distintas diferencias composicionales. La proyección PCA tridimensional (Figura 1F) valida aún más la naturaleza multidimensional de la heterogeneidad celular, y el color de PC3 revela una estratificación adicional más allá de los componentes de varianza primarios.

Caracterización de grupos y análisis de trayectoria funcional

La delimitación detallada de los límites de los grupos mediante el análisis de la envolvente convexa (Figura 2A) establece una clara demarcación entre los subconjuntos de células T funcionalmente distintos. El análisis de velocidad del ARN (Figura 2B) revela estados transcripcionales dinámicos y transiciones direccionales del destino celular, con líneas de flujo de velocidad que indican las principales trayectorias de diferenciación hacia fenotipos de agotamiento. El análisis de la varianza explicada (Figura 2C) demuestra que los tres primeros componentes principales capturan aproximadamente el 10% de la varianza acumulada, con un punto de inflexión en PC8 que sugiere una dimensionalidad óptima. El análisis del gráfico dirigido por fuerza (Figura 2D) ilustra las relaciones interconectadas entre los tipos de células, lo que revela una organización jerárquica y vías de transición. El mapeo de difusión (Figura 2E) captura el panorama de la función citotóxica, identificando regiones de alta y baja citotoxicidad a lo largo de un gradiente continuo. El agrupamiento jerárquico (Figura 2F) categoriza los tipos de células en función de la similitud transcripcional, con un umbral de aproximadamente el 60% de distancia que revela las principales divisiones funcionales en todo el conjunto de datos.

Firmas moleculares y caracterización funcional en los subconjuntos de células T

El perfil de expresión génica en los estados de agotamiento y efector (Figura 3A) revela patrones de agrupamiento jerárquico que distinguen múltiples subtipos de células T CD27+, incluidos los subtipos de agotamiento terminal, efector y similares a la memoria. El mapa de calor demuestra la expresión diferencial de los principales marcadores de agotamiento (PDCD1, CTLA4, HAVCR2) y las moléculas efectoras en los grupos. El análisis de la expresión de CD27 (Figura 3B) muestra niveles significativamente elevados en el tejido tumoral en comparación con el tejido adyacente normal en múltiples subconjuntos de células T, lo que valida la relevancia de las poblaciones CD27+ en la patogénesis del CRC. La matriz de diagrama de puntos (Figura 3C) perfila sistemáticamente la expresión de marcadores específicos de cada tipo de célula, lo que revela firmas moleculares distintas para cada grupo. Los diagramas de violín apilados (Figura 3D) demuestran la distribución de la expresión de CD27 en los subconjuntos de células T CD8+, con una variación notable entre los compartimentos de agotamiento, efector y memoria. El análisis de enriquecimiento de vías (Figura 3E) identifica los principales procesos biológicos asociados con cada subconjunto, incluido la señalización de las citocinas, la reprogramación metabólica y la regulación del ciclo celular. Los diagramas de cresta (Figura 3F) ilustran las distribuciones de expresión de los genes críticos (CD27, CD8A, PRF1, CD3) en los diferentes estados funcionales.

Distribución tisular y elaboración de perfiles del estado funcional

El análisis de clústeres (Figura 4A) revela patrones de firma celular en muestras individuales, lo que demuestra la heterogeneidad específica del paciente en la composición de las células T. El gráfico de barras apiladas (Figura 4B) cuantifica las proporciones de los tipos de células en el tejido adyacente normal frente al tumor, lo que muestra un enriquecimiento drástico de las poblaciones de agotamiento en el microambiente tumoral, con las células efectoras que comprenden la fracción más grande. El análisis del gráfico de radar (Figura 4C) compara los estados funcionales de los subtipos de células T CD8+ en cinco dimensiones clave: agotamiento, citotoxicidad, activación, proliferación y memoria. Las células T CD27+/CD8+ de agotamiento terminal exhiben las puntuaciones de agotamiento más altas, mientras que las células T CD27-/CD8+ efectoras muestran la máxima citotoxicidad. El análisis de la expresión génica diferencial (Figura 4D) identifica los principales genes sobreexpresados y subexpresados, con IGHG1-4 que muestra una sobreexpresión drástica (log2FC > 2). Los diagramas de caja (Figura 4E) demuestran la dinámica de la expresión de los marcadores de agotamiento (GZMB, PDCD1) y las moléculas reguladoras (TOX, TCF7) en los subtipos. La comparación de la puntuación de citotoxicidad (Figura 4F) revela diferencias significativas entre el tejido tumoral y el tejido normal en todos los subconjuntos de CD27+/CD8+, con las poblaciones de agotamiento que muestran una capacidad citotóxica reducida.

Trayectorias de diferenciación y expresión génica dinámica

El análisis del diagrama de volcanes (Figura 5A) destaca los principales genes expresados de forma diferencial entre las poblaciones de agotamiento y efector, con PDCD1, HAVCR2 y LAG3 que muestran una regulación al alza significativa en las células de agotamiento. El diagrama MA (Figura 5B) demuestra la expresión diferencial dependiente de la expresión, lo que identifica tanto los genes de alta abundancia como los de baja abundancia con cambios significativos. El análisis de la trayectoria de pseudotiempo (Figura 5C) revela la progresión del agotamiento de las células T desde los estados de agotamiento terminal hasta los estados de efector, con las células que pasan por fases de agotamiento intermedias. El análisis de la trayectoria de ramificación (Figura 5D) identifica tres destinos de diferenciación distintos: vías de efector, memoria y agotamiento, con puntos de ramificación que indican nodos de decisión críticos en la determinación del destino de las células T. La dinámica de la expresión génica a lo largo de la trayectoria de agotamiento (Figura 5E) muestra patrones temporales distintos para los principales marcadores: CD69, PRF1 y los genes HLA alcanzan su punto máximo en las etapas iniciales, mientras que PDCD1, LAG3 y HAVCR2 aumentan progresivamente hacia el agotamiento terminal. El análisis del campo de velocidad del ARN (Figura 5F) ilustra la dinámica transcripcional direccional subyacente al proceso de agotamiento, con líneas de flujo que convergen hacia los estados de agotamiento terminal.

Transiciones de estado y dinámica de la regulación génica

El análisis de conectividad PAGA (Figura 6A) construye una red integral de transición de estado, lo que revela relaciones jerárquicas entre los estados ingenuo, activado, de memoria, efector y de agotamiento. El gráfico demuestra que las poblaciones de agotamiento representan puntos finales de diferenciación terminal con una plasticidad limitada. El diagrama de flujo de transición del estado celular (Figura 6B) cuantifica las probabilidades de transición entre los estados funcionales, lo que muestra las principales vías desde el estado activado (p = 0,79) a través de los estados preagotamiento y de memoria hasta el agotamiento terminal. El árbol de trayectoria de estilo Monocle (Figura 6C) proporciona una visualización de trayectoria alternativa, lo que confirma múltiples eventos de ramificación que conducen a diversos puntos finales. Las curvas principales de Slingshot (Figura 6D) capturan trayectorias suaves a lo largo de la línea a través del espacio de los componentes principales, con distintas vías para cada ruta de diferenciación. El mapa de calor de la cascada de expresión génica (Figura 6E) revela patrones de regulación temporal a lo largo de la trayectoria de agotamiento, con los genes efectores iniciales (CD69, PRF1) que pasan a la expresión intermedia de puntos de control (PDCD1, CTLA4) y los marcadores de agotamiento terminal (HAVCR2, TOX, LAG3). El análisis de correlación (Figura 6F) demuestra relaciones inversas entre las puntuaciones de agotamiento y la citotoxicidad a lo largo del pseudotiempo (r = -0,905), lo que confirma la disminución funcional durante la progresión del agotamiento.

Redes de regulación génica e interacciones inmunitarias

El análisis GSEA (Figura 7A) revela firmas de agotamiento de células T enriquecidas (NES = 0,8, p = 0,047) con curvas de clasificación que identifican los principales genes contribuyentes, incluidos los impulsores de las células T de memoria y los reguladores del ciclo celular. El análisis de la red de regulación génica (Figura 7B) construye interacciones entre los genes diana de los factores de transcripción que gobiernan el agotamiento de las células T, lo que destaca los principales nodos reguladores, incluidos los factores de transcripción que controlan las transiciones de efector a agotamiento. El mapeo de la interacción ligando-receptor (Figura 7C) en el microambiente tumoral del CRC identifica las principales vías de comunicación entre células, incluido PD-1/PD-L1, CTLA-4/CD80/CD86 y otras interacciones de puntos de control. Los patrones de expresión de citocinas (Figura 7D) revelan una visualización de anillos concéntricos de la producción de citocinas en los tipos de células, con IFN-γ, TNF-α, IL-2 y granzimas que muestran una expresión diferencial. El análisis de la expresión de puntos de control inmunitarios (Figura 7E) demuestra distintos patrones en los subtipos de células T CD8+ para PDCD1, CTLA4, LAG3 y ENTPD1, con diagramas de caja que revelan la mayor expresión en las poblaciones de agotamiento. El análisis composicional (Figura 7F) muestra la distribución de los principales clonotipos de TCR y su diversidad en las muestras, con un gráfico de barras que indica los patrones de expansión clonal y un gráfico circular que revela la composición general de los clonotipos.

Relevancia clínica e implicaciones terapéuticas

La simulación de la distribución espacial (Figura 8A) modela la organización celular dentro de la arquitectura del tejido del CRC, lo que revela que las células T CD27+ de agotamiento se concentran en las regiones centrales del tumor, mientras que las poblaciones de efector se distribuyen más periféricamente. La red de comunicación entre células (Figura 8B) mapea la señalización intercelular entre las células T, las células mieloides y otros componentes inmunitarios, lo que destaca las interacciones de las células T reguladoras y la supresión mediada por las células mieloides. El análisis de supervivencia (Figura 8C) presenta un análisis de supervivencia de Kaplan-Meier exploratorio estratificado por las puntuaciones de agotamiento de las células T. El análisis arrojó HR = 0,505 con p = 0,673, lo que no alcanza la significación estadística. Este hallazgo se presenta explícitamente como una hipótesis generadora; es probable que la cohorte disponible sea de tamaño insuficiente, como lo apoya el cálculo de la potencia posterior. No se demostró ninguna asociación pronóstica estadísticamente significativa, y es necesaria una validación prospectiva en cohortes más grandes antes de que se puedan extraer conclusiones clínicas. El análisis del objetivo terapéutico de expresión (Figura 8D) compara las células T de agotamiento y efector, lo que identifica PDCD1 y LAG3 como las moléculas de puntos de control más sobreexpresadas en las poblaciones de agotamiento, lo que sugiere objetivos de inmunoterapia prioritarios. El panel de heterogeneidad funcional de las células T citotóxicas CD27+ (Figura 8E) integra múltiples análisis: un gráfico circular de la composición del tipo de célula que muestra los principales subconjuntos, histogramas de la puntuación de agotamiento que revelan una distribución bimodal, un diagrama de dispersión de la función en el agotamiento que demuestra una correlación negativa y un mapa de calor de la expresión de marcadores que muestra las firmas. El panel completo (Figura 8F) sintetiza los hallazgos en un mapa del panorama de agotamiento, mapas de calor funcional-composicional que muestran patrones específicos del subconjunto, trayectorias fenotípicas que ilustran las transiciones de estado y curvas de supervivencia estratificadas por los marcadores de agotamiento que demuestran el valor pronóstico clínico.

Expresión diferencial de los marcadores de agotamiento en las líneas celulares de CRC

Este gráfico de barras presenta un análisis comparativo de la expresión de marcadores de agotamiento de las células T entre cuatro líneas celulares de cáncer colorrectal tras el cocultivo con PBMC: las líneas celulares MSI-H HCT116 y RKO, y las líneas celulares MSS HCT15 y SW480. Todos los valores están normalizados con respecto a sus respectivas líneas de base de monocultivo (indicadas por la línea discontinua gris en y = 1). La figura demuestra claramente que las líneas celulares MSI-H exhiben una inducción significativamente mayor de los cuatro marcadores de agotamiento en comparación con las líneas celulares MSS. Para PDCD1, HCT116 mostró una sobreexpresión de 5,2 veces y RKO de 4,8 veces, mientras que SW480 mostró un aumento de 2,1 veces y HCT15 de 1,9 veces. LAG3 sigue un patrón similar: HCT116 (4,8 veces), RKO (4,2 veces), SW480 (1,8 veces) y HCT15 (1,6 veces). HAVCR2 y CD27 muestran diferencias progresivamente menores, pero aún notables, entre las líneas celulares MSI-H y MSS. Los modelos MSI-H (HCT116 y RKO) demuestran consistentemente una mayor respuesta (2-3 veces) al cocultivo con células inmunitarias en todos los marcadores de agotamiento en comparación con los modelos MSS (HCT15 y SW480), lo que apoya la hipótesis de que los cánceres colorrectales inmunogénicos crean un microambiente más permisivo para el desarrollo del fenotipo de agotamiento de las células T (Figura 9).

Arquitectura celular y análisis de reducción de dimensionalidad de las poblaciones de células T CD27+

La caracterización inicial revela patrones de agrupamiento distintos de las células T citotóxicas CD27+ dentro del microambiente tumoral del cáncer colorrectal. El análisis UMAP (Figura 1A) identifica múltiples grupos celulares discretos con diversos estados fenotípicos, lo que demuestra una heterogeneidad sustancial dentro del compartimento de células T CD27+. La visualización t-SNE (Figura 1B) captura un gradiente de agotamiento, con células que pasan de estados funcionales a estados de agotamiento terminal, como lo indica el gradiente de color que representa las puntuaciones de agotamiento. La Figura 1C presenta un análisis de contribución de genes basado en PCA, que revela los principales impulsores moleculares de la heterogeneidad de las células T. La distribución de la densidad celular (Figura 1D) destaca los patrones de enriquecimiento en todo el microambiente tisular, mientras que el análisis comparativo entre el tejido tumoral y el tejido adyacente normal (Figura 1E) demuestra distintas diferencias composicionales. La proyección PCA tridimensional (Figura 1F) valida aún más la naturaleza multidimensional de la heterogeneidad celular, con el color de PC3 que revela una estratificación adicional más allá de los componentes de varianza primarios.

Caracterización de grupos y análisis de trayectoria funcional

La delineación detallada de los límites de los grupos mediante el análisis de la envolvente convexa (Figura 2A) establece una clara demarcación entre los subconjuntos de células T funcionalmente distintos. El análisis de la velocidad del ARN (Figura 2B) revela estados transcripcionales dinámicos y transiciones direccionales del destino celular, con líneas de flujo de velocidad que indican las principales trayectorias de diferenciación hacia fenotipos de agotamiento. El análisis de la varianza explicada (Figura 2C) demuestra que los tres primeros componentes principales capturan aproximadamente el 10% de la varianza acumulada, con un punto de inflexión en PC8 que sugiere una dimensionalidad óptima. El análisis del gráfico dirigido por fuerza (Figura 2D) ilustra las relaciones interconectadas entre los tipos de células, revelando una organización jerárquica y vías de transición. El mapeo de difusión (Figura 2E) captura el paisaje de la función citotóxica, identificando regiones de alta y baja citotoxicidad a lo largo de un gradiente continuo. El agrupamiento jerárquico (Figura 2F) categoriza los tipos de células en función de la similitud transcripcional, con un umbral de aproximadamente el 60% de distancia que revela las principales divisiones funcionales en todo el conjunto de datos.

Firmas moleculares y caracterización funcional en los subconjuntos de células T

El perfil de expresión génica en los estados de agotamiento y efector (Figura 3A) revela patrones de agrupamiento jerárquico que distinguen múltiples subtipos de células T CD27+, incluidos los subtipos terminales de agotamiento, efectores y similares a la memoria. El mapa de calor demuestra la expresión diferencial de los principales marcadores de agotamiento (PDCD1, CTLA4, HAVCR2) y las moléculas efectoras en los grupos. El análisis de la expresión de CD27 (Figura 3B) muestra niveles significativamente elevados en el tejido tumoral en comparación con el tejido adyacente normal en múltiples subconjuntos de células T, lo que valida la relevancia de las poblaciones CD27+ en la patogénesis del CRC. La matriz de diagrama de puntos (Figura 3C) perfila sistemáticamente la expresión de marcadores específicos de cada tipo de célula, revelando firmas moleculares distintas para cada grupo. Los diagramas de violín apilados (Figura 3D) demuestran la distribución de la expresión de CD27 en los subconjuntos de células T CD8+, con una variación notable entre los compartimentos de agotamiento, efector y memoria. El análisis de enriquecimiento de vías (Figura 3E) identifica los principales procesos biológicos asociados con cada subconjunto, incluidos la señalización de citocinas, la reprogramación metabólica y la regulación del ciclo celular. Los diagramas de cresta (Figura 3F) ilustran las distribuciones de expresión de genes críticos (CD27, CD8A, PRF1, CD3) en diferentes estados funcionales.

Distribución tisular y perfilado del estado funcional

El análisis de Clustermap (Figura 4A) revela patrones de firma celular en muestras individuales, lo que demuestra la heterogeneidad específica del paciente en la composición de las células T. El diagrama de barras apilado (Figura 4B) cuantifica las proporciones de tipos de células en el tejido adyacente normal frente al tumor, lo que muestra un enriquecimiento drástico de las poblaciones de agotamiento en el microambiente tumoral, con las células efectoras que comprenden la fracción más grande. El análisis del diagrama de radar (Figura 4C) compara los estados funcionales de los subtipos de células T CD8+ en cinco dimensiones clave: agotamiento, citotoxicidad, activación, proliferación y memoria. Las células T CD27+/CD8+ de agotamiento terminal exhiben las puntuaciones de agotamiento más altas, mientras que las células T CD27-/CD8+ efectoras muestran la máxima citotoxicidad. El análisis de la expresión génica diferencial (Figura 4D) identifica los principales genes sobreexpresados y subexpresados, con IGHG1-4 que muestra una sobreexpresión drástica (log2FC > 2). Los diagramas de caja (Figura 4E) demuestran la dinámica de la expresión de los marcadores de agotamiento (GZMB, PDCD1) y las moléculas reguladoras (TOX, TCF7) en los subtipos. La comparación de la puntuación de citotoxicidad (Figura 4F) revela diferencias significativas entre el tejido tumoral y el tejido adyacente normal en todos los subconjuntos de CD27+/CD8+, con las poblaciones de agotamiento que muestran una capacidad citotóxica reducida.

Trayectorias de diferenciación y expresión génica dinámica

El análisis del diagrama de volcanes (Figura 5A) destaca los principales genes expresados de forma diferencial entre las poblaciones de agotamiento y efector, con PDCD1, HAVCR2 y LAG3 que muestran una regulación al alza significativa en las células de agotamiento. El diagrama MA (Figura 5B) demuestra la expresión dependiente de la expresión diferencial, identificando tanto los genes de alta abundancia como los de baja abundancia con cambios significativos. El análisis de la trayectoria de pseudotiempo (Figura 5C) revela la progresión del agotamiento de las células T desde los estados de agotamiento terminal hasta los estados de efector, con las células que pasan por fases de agotamiento intermedias. El análisis de la trayectoria de ramificación (Figura 5D) identifica tres destinos de diferenciación distintos: vías de efector, memoria y agotamiento, con puntos de ramificación que indican nodos de decisión críticos en la determinación del destino de las células T. La dinámica de la expresión génica a lo largo de la trayectoria de agotamiento (Figura 5E) muestra patrones temporales distintos para los principales marcadores: CD69, PRF1, los genes HLA alcanzan su punto máximo en las primeras etapas, mientras que PDCD1, LAG3 y HAVCR2 aumentan progresivamente hacia el agotamiento terminal. El análisis del campo de velocidad del ARN (Figura 5F) ilustra la dinámica transcripcional direccional subyacente al proceso de agotamiento, con líneas de flujo que convergen hacia los estados de agotamiento terminal.

Transiciones de estado y dinámica de la regulación génica

El análisis de conectividad PAGA (Figura 6A) construye una red integral de transición de estado, revelando relaciones jerárquicas entre los estados ingenuo, activado, de memoria, efector y de agotamiento. El gráfico demuestra que las poblaciones de agotamiento representan puntos finales de diferenciación terminal con una plasticidad limitada. El diagrama de flujo de transición del estado celular (Figura 6B) cuantifica las probabilidades de transición entre los estados funcionales, mostrando las principales vías desde el estado activado (p = 0,79) a través de los estados pre-agotamiento y de memoria hacia el agotamiento terminal. El árbol de trayectoria de estilo Monocle (Figura 6C) proporciona una visualización de trayectoria alternativa, confirmando múltiples eventos de ramificación que conducen a diversos puntos finales. Las curvas principales de Slingshot (Figura 6D) capturan trayectorias suaves a lo largo de la línea a través del espacio de componentes principales, con distintas vías para cada ruta de diferenciación. El mapa de calor de la cascada de expresión génica (Figura 6E) revela patrones de regulación temporal a lo largo de la trayectoria de agotamiento, con los primeros genes de efector (CD69, PRF1) que pasan a la expresión intermedia de puntos de control (PDCD1, CTLA4) y los marcadores de agotamiento terminal (HAVCR2, TOX, LAG3). El análisis de correlación (Figura 6F) demuestra relaciones inversas entre las puntuaciones de agotamiento y la citotoxicidad a lo largo del pseudotiempo (r = -0,905), lo que confirma la disminución funcional durante la progresión del agotamiento.

Redes de regulación génica e interacciones inmunitarias

El análisis GSEA (Figura 7A) revela firmas de agotamiento de células T enriquecidas (NES = 0,8, p = 0,047) con curvas de clasificación que identifican los principales genes contribuyentes, incluidos los impulsores de las células T de memoria y los reguladores del ciclo celular. El análisis de la red de regulación génica (Figura 7B) construye interacciones entre los factores de transcripción y los genes diana que gobiernan el agotamiento de las células T, destacando los principales nodos reguladores, incluidos los factores de transcripción que controlan las transiciones de efector a agotamiento. El mapeo de la interacción ligando-receptor (Figura 7C) en el microambiente tumoral del CRC identifica las principales vías de comunicación entre células, incluido PD-1/PD-L1, CTLA-4/CD80/CD86 y otras interacciones de puntos de control. Los patrones de expresión de citocinas (Figura 7D) revelan una visualización de anillos concéntricos de la producción de citocinas en los tipos de células, con IFN-γ, TNF-α, IL-2 y granzimas que muestran una expresión diferencial. El análisis de la expresión de puntos de control inmunitarios (Figura 7E) demuestra distintos patrones en los subtipos de células T CD8+ para PDCD1, CTLA4, LAG3 y ENTPD1, con diagramas de caja que revelan la mayor expresión en las poblaciones de agotamiento. El análisis composicional (Figura 7F) muestra la distribución de los principales clonotipos de TCR y su diversidad en las muestras, con un diagrama de barras que indica los patrones de expansión clonal y un diagrama circular que revela la composición general de los clonotipos.

Relevancia clínica e implicaciones terapéuticas

Simulación de la distribución espacial (Figura 8A) que modela la organización celular dentro de la arquitectura del tejido del CCRC, revelando que las células T CD27+ exhaustas se concentran en las regiones del núcleo tumoral, mientras que las poblaciones efectoras se distribuyen más periféricamente. El mapa de la red de comunicación célula-célula (Figura 8B) muestra la señalización intercelular entre las células T, las células mieloides y otros componentes inmunitarios, destacando las interacciones de las células T reguladoras y la supresión mediada por las células mieloides. El análisis de supervivencia (Figura 8C) presenta un análisis exploratorio de supervivencia de Kaplan-Meier estratificado por las puntuaciones de agotamiento de las células T. El análisis arrojó un HR = 0,505 con p = 0,673, lo que no alcanza significación estadística. Este hallazgo se presenta explícitamente como una hipótesis que debe ser investigada; la cohorte disponible probablemente es demasiado pequeña, como lo demuestra el cálculo de la potencia post-hoc. No se demostró ninguna asociación pronóstica estadísticamente significativa, y es necesaria una validación prospectiva en cohortes más grandes antes de que se puedan extraer conclusiones clínicas. El análisis de la expresión de los objetivos terapéuticos (Figura 8D) compara las células T exhaustas con las células T efectoras, identificando PDCD1 y LAG3 como las moléculas de control inmunitario más significativamente sobreexpresadas en las poblaciones exhaustas, lo que sugiere objetivos de inmunoterapia prioritarios. El panel de heterogeneidad funcional de las células T citotóxicas CD27+ (Figura 8E) integra múltiples análisis: un gráfico circular de la composición del tipo de célula que muestra los subconjuntos dominantes, histogramas de la puntuación de agotamiento que revelan una distribución bimodal, un diagrama de dispersión de la función en relación con el agotamiento que demuestra una correlación negativa y un mapa de calor de la expresión génica clave que muestra las firmas de marcadores. El panel completo (Figura 8F) sintetiza los hallazgos en un mapa del paisaje de agotamiento, mapas de calor funcional-compositivos que muestran patrones específicos de subconjuntos, trayectorias fenotípicas que ilustran las transiciones de estado y curvas de supervivencia estratificadas por marcadores de agotamiento que demuestran el valor pronóstico clínico.

Expresión diferencial de los marcadores de agotamiento en las líneas celulares de CCRC

Este gráfico de barras presenta un análisis comparativo de la expresión de los marcadores de agotamiento de las células T entre cuatro líneas celulares de cáncer colorrectal después del cocultivo con PBMCs: las líneas celulares MSI-H HCT116 y RKO, y las líneas celulares MSS HCT15 y SW480. Todos los valores están normalizados con respecto a sus respectivas líneas de base de monocultivo (indicado por la línea discontinua gris en y = 1). La figura demuestra claramente que las líneas celulares MSI-H exhiben una inducción sustancialmente mayor de los cuatro marcadores de agotamiento en comparación con las líneas celulares MSS. Para PDCD1, HCT116 mostró una sobreexpresión de 5,2 veces y RKO mostró una sobreexpresión de 4,8 veces, mientras que SW480 mostró un aumento de 2,1 veces y HCT15 mostró un aumento de 1,9 veces. LAG3 sigue un patrón similar: HCT116 (4,8 veces), RKO (4,2 veces), SW480 (1,8 veces) y HCT15 (1,6 veces). HAVCR2 y CD27 muestran diferencias progresivamente menores, pero aún notables, entre las líneas celulares MSI-H y MSS. Los modelos MSI-H (HCT116 y RKO) demuestran consistentemente una mayor respuesta de 2 a 3 veces al cocultivo con células inmunitarias en todos los marcadores de agotamiento en comparación con los modelos MSS (HCT15 y SW480), lo que apoya la hipótesis de que los cánceres colorrectales inmunogénicos crean un microambiente más permisivo para el desarrollo del fenotipo de agotamiento de las células T (Figura 9).

Discusión

En este estudio transcriptómico unicelular integral, hemos dilucidado la heterogeneidad celular, los estados funcionales y la dinámica de agotamiento de las células T citotóxicas CD27+ dentro del microambiente tumoral del cáncer colorrectal. Nuestros hallazgos revelan 11 grupos distintos de células T CD27+ con diversas características fenotípicas, trayectorias de agotamiento progresivas y patrones de organización espacial que proporcionan información crítica sobre los mecanismos subyacentes a la disfunción inmunitaria en el CCRC e identifican posibles vulnerabilidades terapéuticas para la optimización de la inmunoterapia.

Nuestros análisis UMAP y t-SNE identificaron una heterogeneidad sustancial dentro del compartimento de células T citotóxicas CD27+, que abarca poblaciones terminales exhaustas, efectoras y similares a la memoria con distintas firmas moleculares y capacidades funcionales. Este hallazgo se alinea con estudios unicelulares recientes que demuestran que los linfocitos infiltrantes tumorales existen a lo largo de un continuo de estados funcionales en lugar de categorías discretas ([Guo et al., 2018]). La identificación de múltiples subconjuntos exhaustos dentro de las células T CD27+ es particularmente significativa, ya que la evidencia reciente sugiere que las células T exhaustas comprenden poblaciones heterogéneas con diferente capacidad de respuesta al bloqueo de los puntos de control inmunitario ([Miller et al., 2019]). Específicamente, las células T progenitoras exhaustas que expresan TCF1 conservan la capacidad proliferativa y pueden ser revitalizadas por el bloqueo de PD-1, mientras que las células T terminalmente exhaustas exhiben una disfunción irreversible ([Im et al., 2016]). Nuestro análisis de trayectoria reveló que las células T CD27+ progresan a través de fases intermedias de agotamiento antes de alcanzar estados de agotamiento terminal, lo que sugiere posibles puntos de intervención para prevenir la disfunción irreversible.

La expresión elevada de CD27 en el tejido tumoral en comparación con el tejido adyacente normal en múltiples subconjuntos de células T valida la relevancia biológica de esta población en la patogénesis del CCRC. Se ha demostrado que la señalización de CD27 promueve la supervivencia de las células T y mejora las funciones efectoras a través de la activación de NF-κB y la reprogramación metabólica ([Borst et al., 2005]). Sin embargo, la presencia paradójica de células CD27+ dentro de las poblaciones exhaustas sugiere una dinámica regulatoria compleja dentro del microambiente tumoral (TME). La base mecánica de esta paradoja se comprende cada vez más: en condiciones de coestimulación aguda, la activación de NF-κB mediada por CD27 promueve programas de genes pro-supervivencia y efectores. Sin embargo, en la exposición crónica a antígenos en el TME, la señalización sostenida de CD27 impulsa la actividad persistente de NF-κB a través de bucles de transcripción de retroalimentación, aumentando paradójicamente los programas de genes pro-apoptóticos y pro-agotamiento, incluidos BLIMP-1 y los miembros de la familia NR4A. Al mismo tiempo, la señalización coestimuladora prolongada se cruza con una disfunción mitocondrial progresiva, caracterizada por la pérdida de la capacidad respiratoria de reserva, el aumento de la producción de especies reactivas de oxígeno (ROS) y la alteración de la fosforilación oxidativa, alteraciones metabólicas que se han relacionado mecánicamente con la progresión del agotamiento en modelos de estimulación antigénica crónica. Estas alteraciones metabólicas y de señalización ascendentes convergen en última instancia en la remodelación de la cromatina mediada por TOX, consolidando epigenéticamente el estado de agotamiento de manera en gran medida irreversible. Estudios recientes han demostrado que la señalización crónica de CD27 en presencia de una estimulación antigénica persistente puede contribuir paradójicamente al agotamiento de las células T a través del estrés metabólico sostenido y la disfunción mitocondrial ([Peperzak et al., 2010]). Este doble papel de CD27 tanto en la promoción de la función de las células T como en la posible contribución al agotamiento en condiciones de estimulación crónica justifica una mayor investigación y puede tener implicaciones importantes para las estrategias terapéuticas dirigidas a CD27.

Nuestra caracterización molecular integral identificó una regulación progresiva de múltiples receptores inhibitorios, incluidos PDCD1 (PD-1), LAG3, HAVCR2 (TIM-3) y CTLA4 a lo largo de la trayectoria de agotamiento, acompañada de una disminución de las moléculas efectoras, como PRF1 y GZMB. Esta coexpresión de múltiples moléculas de control inmunitario, denominada "inhibición en capas", representa una característica fundamental de las células T severamente exhaustas y contribuye a una alteración funcional profunda ([Blackburn et al., 2008]). La fuerte correlación inversa (r = -0,905) entre las puntuaciones de agotamiento y la citotoxicidad a lo largo de la progresión del pseudotiempo demuestra cuantitativamente las consecuencias funcionales del desarrollo del agotamiento. Es importante destacar que nuestro análisis de la cascada de expresión génica reveló distintas fases temporales de agotamiento, con los primeros genes efectores (CD69, PRF1) que hacen la transición a la expresión intermedia de puntos de control (PDCD1, CTLA4) y, en última instancia, a los marcadores de agotamiento terminal (HAVCR2, TOX, LAG3). Este patrón temporal sugiere que el agotamiento se desarrolla a través de programas de transcripción coordinados en lugar de cambios en la expresión génica estocástica.

El análisis de la red reguladora de genes identificó los principales factores de transcripción que gobiernan las transiciones de efector a exhausto, destacando el papel central de la reprogramación transcripcional en el desarrollo del agotamiento. TOX ha surgido como un regulador maestro del agotamiento de las células T, con estudios que demuestran que TOX es tanto necesario como suficiente para establecer el programa de agotamiento a través de la remodelación de la cromatina y las modificaciones epigenéticas ([Scott et al., 2019]). Nuestro hallazgo de una regulación progresiva de TOX a lo largo de la trayectoria de agotamiento apoya su papel fundamental en la disfunción de las células T CD27+ en el CCRC. Además, la compleja interacción entre TCF1 y TOX determina el equilibrio entre los estados de agotamiento progenitor y terminal, con las células progenitoras TCF1+ que mantienen la capacidad de autorrenovación, mientras que las células terminalmente exhaustas TOXhi exhiben una disfunción fija ([Siddiqui et al., 2019]). Comprender estas jerarquías transcripcionales proporciona información mecánica sobre el establecimiento del agotamiento y las posibles estrategias para la reprogramación epigenética de las células T exhaustas.

Nuestros experimentos de cocultivo in vitro utilizando cuatro líneas celulares de cáncer colorrectal demostraron que las líneas celulares MSI-H (HCT116 y RKO) indujeron de 2 a 3 veces más expresión de marcadores de agotamiento en comparación con las líneas celulares MSS (HCT15 y SW480). Se observaron diferencias particularmente dramáticas en la inducción de PDCD1 (HCT116: 5,2 veces, RKO: 4,8 veces frente a SW480: 2,1 veces, HCT15: 1,9 veces) y la inducción de LAG3 (HCT116: 4,8 veces, RKO: 4,2 veces frente a SW480: 1,8 veces, HCT15: 1,6 veces). Estos hallazgos se basan en la cuantificación de ARNm mediante qRT-PCR y deben interpretarse como evidencia a nivel de la transcripción consistente con, en lugar de confirmar completamente, las firmas de agotamiento computacionales. Los mecanismos reguladores postranscripcionales pueden desacoplar la abundancia de transcritos y proteínas; la confirmación a nivel de proteínas mediante el tinte de superficie por citometría de flujo para PD-1 y LAG-3 representa el siguiente paso experimental lógico y se identifica como una prioridad para el trabajo futuro. A pesar de esta limitación, el patrón diferencial entre las líneas celulares MSI-H y MSS es sólido en cuatro marcadores independientes y tres réplicas biológicas, lo que proporciona una consistencia interna a los hallazgos transcripcionales. Este hallazgo proporciona evidencia experimental directa de la capacidad inmunogénica diferencial entre los tumores MSI-H y MSS y ayuda a explicar las diferentes respuestas clínicas al bloqueo de los puntos de control inmunitario observadas en los pacientes con CCRC. Los tumores MSI-H albergan altas cargas de mutación que resultan de un mecanismo de reparación de errores defectuoso, lo que conduce a un aumento de la producción de neoantígenos y una mayor infiltración de células T ([Boland y Goel, 2010]). La elevada carga de neoantígenos en los tumores MSI-H impulsa la activación crónica de las células T y el agotamiento posterior, creando paradójicamente tanto la justificación como el desafío para la inmunoterapia ([Germano et al., 2017]).

La observación de que los tumores MSI-H inducen una mayor expresión de los marcadores de agotamiento sugiere que estos tumores crean un microambiente inmunológico más "caliente" con un compromiso inmunitario activo, aunque acompañado de una regulación ascendente compensatoria de las vías inhibitorias. Esto contrasta con los tumores MSS, que exhiben fenotipos excluidos por el sistema inmunitario o desérticos, caracterizados por una escasa infiltración de células T y una activación inmunitaria limitada ([Becht et al., 2016]). Es importante destacar que, si bien los tumores MSI-H muestran un mayor agotamiento de base, también demuestran respuestas superiores al bloqueo de PD-1/PD-L1, lo que sugiere que revertir el agotamiento en un sistema inmunitario ya comprometido es más factible que superar la exclusión inmunitaria fundamental ([Overman et al., 2017]). Estos hallazgos subrayan la importancia de considerar los subtipos moleculares de los tumores al diseñar estrategias inmunoterapéuticas y sugieren que los tumores MSS pueden requerir enfoques combinados para mejorar simultáneamente la infiltración de células T y prevenir el agotamiento.

Nuestro análisis computacional de la distribución espacial reveló que las células T CD27+ exhaustas se concentran en las regiones centrales del tumor, mientras que las poblaciones efectoras se distribuyen más periféricamente, lo que destaca la heterogeneidad espacial de las respuestas inmunitarias dentro de los tejidos de CRC. Esta compartimentación espacial refleja el gradiente inmunosupresor dentro de los tumores, donde las regiones centrales exhiben concentraciones más altas de señales inhibitorias, competencia metabólica e hipoxia, que en conjunto promueven la disfunción de las células T ([Salmon et al., 2012]). Estudios recientes de transcriptómica espacial han demostrado que el agotamiento de las células T no se distribuye uniformemente en los tejidos tumorales, sino que se correlaciona con la proximidad a las células malignas, lo que sugiere que los factores microambientales locales regulan de manera crítica los estados funcionales de las células T ([Jerby-Arnon et al., 2018]).

Limitaciones del estudio

Varias limitaciones de este estudio merecen consideración. En primer lugar, el análisis de secuenciación de ARN de una sola célula proporciona instantáneas transcriptómicas, pero no captura los cambios temporales dinámicos ni las respuestas funcionales a las perturbaciones. Los estudios futuros que empleen un muestreo longitudinal o el seguimiento in vivo de las poblaciones de células T proporcionarían información valiosa sobre la cinética del agotamiento y las respuestas terapéuticas. En segundo lugar, si bien nuestros experimentos de cocultivo in vitro validaron la inducción diferencial de marcadores de agotamiento entre las líneas celulares MSI-H (HCT116 y RKO) y las líneas celulares MSS (HCT15 y SW480), los sistemas de cultivo bidimensionales simplificados no pueden recapitular completamente la compleja arquitectura tridimensional y las interacciones multicelulares del microambiente tumoral. Los sistemas de organoides avanzados o los modelos de ratón humanizados proporcionarían plataformas de validación más relevantes desde el punto de vista fisiológico. En tercer lugar, este estudio se basa exclusivamente en un único conjunto de datos disponible públicamente (GSE144735) para todos los análisis computacionales. Si bien GSE144735 abarca tejidos tumorales y tejidos adyacentes normales correspondientes de múltiples pacientes con CRC, lo que proporciona heterogeneidad interna entre las muestras, la ausencia de una cohorte de validación externa independiente limita sustancialmente la generalización de los hallazgos computacionales. La validación cruzada externa en un conjunto de datos independiente de secuenciación de ARN de una sola célula de CRC se identifica como una prioridad para el trabajo futuro. En cuarto lugar, toda la validación experimental de la inducción de marcadores de agotamiento se realizó exclusivamente a nivel de ARNm mediante qRT-PCR. Dado que los mecanismos reguladores postranscripcionales y traduccionales pueden desacoplar la abundancia de transcritos de la expresión de proteínas, la confirmación a nivel de proteínas de la regulación al alza de PDCD1 (PD-1) y LAG3 mediante tinción de superficie por citometría de flujo o Western blot proporcionaría una capa adicional de evidencia crítica; esto se identifica como el siguiente paso experimental lógico. En quinto lugar, el análisis de distribución espacial presentado en la Figura 8A se deriva completamente de una simulación computacional basada en datos de transcriptómica espacial publicados en CRC, en lugar de un perfil espacial empírico directo. Este modelo espacial simulado no debe interpretarse como evidencia experimental directa; la validación espacial experimental mediante transcriptómica espacial (por ejemplo, 10x Visium) o inmunofluorescencia multiplex en secciones de tejido de CRC sería necesaria para confirmar el patrón espacial inferido. En sexto lugar, el análisis de supervivencia exploratorio (HR = 0,505, p = 0,673) no alcanzó una significación estadística, y el análisis de potencia post hoc indica que la cohorte disponible está sustancialmente infrautilizada para detectar un efecto de esta magnitud. Estos hallazgos se presentan únicamente como generadores de hipótesis y requieren una validación prospectiva en cohortes más grandes y adecuadamente dimensionadas antes de que se puedan extraer conclusiones pronósticas.

Limitaciones del estudio

Varias limitaciones de este estudio merecen consideración. En primer lugar, el análisis de secuenciación de ARN de una sola célula proporciona instantáneas transcriptómicas, pero no captura los cambios temporales dinámicos ni las respuestas funcionales a las perturbaciones. Los estudios futuros que empleen un muestreo longitudinal o el seguimiento in vivo de las poblaciones de células T proporcionarían información valiosa sobre la cinética del agotamiento y las respuestas terapéuticas. En segundo lugar, si bien nuestros experimentos de cocultivo in vitro validaron la inducción diferencial de marcadores de agotamiento entre las líneas celulares MSI-H (HCT116 y RKO) y las líneas celulares MSS (HCT15 y SW480), los sistemas de cultivo bidimensionales simplificados no pueden recapitular completamente la compleja arquitectura tridimensional y las interacciones multicelulares del microambiente tumoral. Los sistemas de organoides avanzados o los modelos de ratón humanizados proporcionarían plataformas de validación más relevantes desde el punto de vista fisiológico. En tercer lugar, este estudio se basa exclusivamente en un único conjunto de datos disponible públicamente (GSE144735) para todos los análisis computacionales. Si bien GSE144735 abarca tejidos tumorales y tejidos adyacentes normales correspondientes de múltiples pacientes con CRC, lo que proporciona heterogeneidad interna entre las muestras, la ausencia de una cohorte de validación externa independiente limita sustancialmente la generalización de los hallazgos computacionales. La validación cruzada externa en un conjunto de datos independiente de secuenciación de ARN de una sola célula de CRC se identifica como una prioridad para el trabajo futuro. En cuarto lugar, toda la validación experimental de la inducción de marcadores de agotamiento se realizó exclusivamente a nivel de ARNm mediante qRT-PCR. Dado que los mecanismos reguladores postranscripcionales y traduccionales pueden desacoplar la abundancia de transcritos de la expresión de proteínas, la confirmación a nivel de proteínas de la regulación al alza de PDCD1 (PD-1) y LAG3 mediante tinción de superficie por citometría de flujo o Western blot proporcionaría una capa adicional de evidencia crítica; esto se identifica como el siguiente paso experimental lógico. En quinto lugar, el análisis de distribución espacial presentado en la Figura 8A se deriva completamente de una simulación computacional basada en datos de transcriptómica espacial publicados en CRC, en lugar de un perfil espacial empírico directo. Este modelo espacial simulado no debe interpretarse como evidencia experimental directa; la validación espacial experimental mediante transcriptómica espacial (por ejemplo, 10x Visium) o inmunofluorescencia multiplex en secciones de tejido de CRC sería necesaria para confirmar el patrón espacial inferido. En sexto lugar, el análisis de supervivencia exploratorio (HR = 0,505, p = 0,673) no alcanzó una significación estadística, y el análisis de potencia post hoc indica que la cohorte disponible está sustancialmente infrautilizada para detectar un efecto de esta magnitud. Estos hallazgos se presentan únicamente como generadores de hipótesis y requieren una validación prospectiva en cohortes más grandes y adecuadamente dimensionadas antes de que se puedan extraer conclusiones pronósticas.

Conclusión

Este estudio transcriptómico integral de una sola célula proporciona un atlas sin precedentes y de alta resolución de la heterogeneidad de las células T citotóxicas CD27+, la dinámica del agotamiento y la organización espacial dentro del microambiente tumoral del cáncer colorrectal. Nuestros hallazgos revelan 11 grupos celulares distintos con trayectorias de diferenciación progresivas que conducen a estados de agotamiento terminales caracterizados por la regulación al alza de múltiples receptores inhibitorios y el deterioro funcional.