Abstract Background Tertiary lymphoid structures (TLS) within the tumor microenvironment have been associated with cancer prognosis and therapeutic response. However, the immunological pattern of a high peritumoral TLS (pTLS) density and its clinical potential in hepatocellular carcinoma (HCC) remain poor. This study aimed to elucidate biological differences related to pTLS density and develop a radiomic classifier for predicting pTLS density in HCC, offering new insights for clinical diagnosis and treatment. Methods Spatial transcriptomics (n=4) and RNA sequencing data (n=952) were used to identify critical regulators of pTLS density and evaluate their prognostic significance in HCC. Baseline MRI images from 660 patients with HCC who had undergone surgery treatment between October 2015 and January 2023 were retrospectively recruited for model development and validation. This included training (n=307) and temporal validation (n=76) cohorts from Xiangya Hospital, and external validation cohorts from three independent hospitals (n=277). Radiomic features were extracted from intratumoral and peritumoral regions of interest and analyzed using machine learning algorithms to develop a predictive classifier. The classifier’s performance was evaluated using the area under the curve (AUC), with prognostic and predictive value assessed across four independent cohorts and in a dual-center outcome cohort of 41 patients who received immunotherapy. Results Patients with HCC and a high pTLS density experienced prolonged median overall survival (p<0.05) and favorable immunotherapy response (p=0.03). Moreover, immune infiltration by mature B cells was observed in the high pTLS density region. Spatial pseudotime analysis and immunohistochemistry staining revealed that expansion of pTLS in HCC was associated with elevated CXCL9 and CXCL10 co-expression. We developed an optimal radiomic-based classifier with excellent discrimination for predicting pTLS density, achieving an AUC of 0.91 (95% CI 0.87, 0.94) in the external validation cohort. This classifier also exhibited promising stratification ability in terms of overall survival (p<0.01), relapse-free survival (p<0.05), and immunotherapy response (p<0.05). Conclusion We identified key regulators of pTLS density in patients with HCC and proposed a non-invasive radiomic classifier capable of assisting in stratification for prognosis and treatment. Keywords: Hepatocellular Carcinoma, Biomarker __________________________________________________________________ WHAT IS ALREADY KNOWN ON THIS TOPIC. WHAT THIS STUDY ADDS * We systematically evaluate the clinical outcomes, spatial characteristics, and formation of high pTLS density in hepatocellular carcinoma. Moreover, we propose a baseline magnetic resonance imaging-based radiomic classifier to predict the high pTLS density. HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OF POLICY * Our study identifies potential novel targets of pTLS density, shedding light on its clinical application in immunotherapy. Developing a non-invasive method for assessing pTLS density could aid in clinical decision-making and treatment selection in hepatocellular carcinoma. Introduction Hepatocellular carcinoma (HCC) is the second leading cause of cancer-related deaths worldwide.[49]^1 Although surgical resection significantly prolongs the median survival time (mST) of patients with HCC, the benefits are often limited by a high recurrence rate of 70% and rapid cancer progression.[50]^2 Nowadays, immune-related molecules, circulating tumor markers and tumor heterogeneity-related gene expression are used as biomarkers for prognosis and treatment response evaluation; however, their accuracy and generalizability require validation through multicenter studies.[51]^1 3 Effective single tumor-specific targets that accurately reflect the biological patterns of HCC are still lacking, prompting interest in spatial structural indicators as emerging biomarkers.[52]^3 4 Tertiary lymphoid structures (TLS) exhibit morphological features and anticancer immune cell infiltration patterns akin to those of lymphoid organs, which are associated with improved clinical outcomes across solid tumors.[53]4,[54]6 In HCC, the presence of intratumoral TLS is associated with favorable prognosis and improved immunotherapy response,[55]7,[56]10 indicating that intratumoral TLS may serve as a promising biomarker to predict HCC outcomes. Similarly, high peritumoral TLS (pTLS) density has demonstrated predictive potential.[57]^11 However, the molecular and spatial characteristics of high pTLS density in HCC require further investigation. Pretherapeutic characterization of high pTLS may facilitate the development of individualized treatment plans for patients with HCC. While the histopathological examination is the gold standard for TLS detection,[58]^4 7 12 a biopsy is invasive and associated with various complications, including bleeding, infection, perforation, and needle tract seeding. Radiological methods offer a non-invasive alternative for disease evaluation and management. Recent studies have demonstrated that radiological features and radiomics could identify the presence of TLS in diverse solid tumors.[59]713,[60]15 Two of these studies investigated the predictive value of CT for intratumoral TLS status. However, compared with CT, MRI provides superior soft-tissue resolution and avoids radiation exposure. Yet, to the best of our knowledge, no studies have applied radiomics to MRI images for predicting pTLS in HCC. Therefore, this study aimed to assess the predictive value of high pTLS density for prognosis and treatment response across multicenter cohorts. Additionally, we investigated the pTLS spatial expression patterns and formation using spatial transcriptomics and validated our findings using immunohistochemical staining. Ultimately, we aimed to develop an MRI-based radiomic model for predicting pTLS to facilitate future non-invasive clinical evaluation. Materials and methods Specimens and clinical data A total of 701 eligible patients with HCC who underwent curative resection were consecutively recruited in four well-known independent hospitals in China from October 2015, to January 2023. Data collected comprised including: (1) Xiangya Hospital data sets, (2) three external validation data sets, and (3) dual-center outcome data sets. The inclusion criterion was histologically confirmed HCC. The exclusion criteria were as follows: (1) poor-quality MRI data (eg, motion artifacts) and pathological slices (eg, absence of tumor site, tissue folds, and fade staining); (2) MRI performed >1 month before surgery; and (3) additional treatment before surgery (eg, liver transplantation, local ablation, or transarterial chemoembolization). Postoperative follow-up involved radiological examination at 3–6-month intervals, tracking clinical outcomes (tumor recurrence, metastasis, and death), with the last follow-up endpoint being June 30, 2023. All clinical features, including age, sex, immediate preoperative alpha-fetoprotein, hepatitis virus infection status (hepatitis B and C virus), albumin, total bilirubin, aspartate transaminase, and alanine transaminase levels, were obtained within 14 days before hepatectomy. Other prognosis-associated variables were collected for multivariable Cox regression analysis, including resection margin and vascular invasion (micro/macro). Additionally, data from 952 patients with HCC were obtained from international multicenter cohorts (The Cancer Genome Atlas-Liver HCC (TCGA-LIHC, n=363; [61]GSE14520, n=242; and [62]GSE76427, n=115); the International Cancer Genome Consortium-Liver Cancer RIKEN (ICGC-LIRI, n=232), and several anti-programmed cell death 1/programmed cell death 1 ligand 1 (anti-PD(L)-1) treatment in solid tumor cohorts (HCC, n=65; non-small cell lung cancer (NSCLC), n=28; pan-cancer, n=100; and urothelial bladder carcinoma (BLCA), n=298). Treatment data In the dual-center outcome data sets, all patients received combination therapy with antiangiogenic agents (lenvatinib 8 or 12 mg; cabozantinib 40 mg; bevacizumab 15 mg/kg every 3 weeks; or ipilimumab 1 mg/kg every 3 weeks) and an intravenously administered anti-PD(L)-1 antibodies[63]^16 (pembrolizumab, 200 mg; atezolizumab, 1,200 mg every 3 weeks, nivolumab, 3 mg/kg, or sintilimab 200 mg). All patients received consistent treatment and were closely monitored. Tumor response was assessed according to the Response Evaluation Criteria in Solid Tumors.[64]^17 Eligible patients were those with a minimum 2-month interval between any previous therapy (eg, locoregional therapy or hepatectomy) and the initiation of combination therapy, and those who had at least one therapeutic response assessment after initiating combination therapy. TLS quantification TLS were identified as ectopic lymphoid structures with lymphoid cells aggregated (especially B and T cells) that lacked integrated organized structures such as capsules.[65]^4 12 18 Independent of clinical and MRI data, pTLS were morphologically detected using H&E-stained images via an automated pipeline ([66]https://github.com/YuMeng-W/TumSeg-main), which achieved a good agreement with manual TLS counting by pathologists using digital H&E-stained images.[67]^19 Briefly, the area of the pre-segmented peritumoral region was calculated, and lymphoid aggregates with a diameter or major axis of the ellipse exceeding 150 µm were extracted.[68]^19 20 The pTLS density was calculated as number/mm^2 in the peritumoral region (defined as 5 mm from the infiltrative tumor border). Patients were categorized as having a high pTLS density (pTLS^high) if the distribution of pTLS was >0.28 number/mm^2, or as having low pTLS density (pTLS^low) if it was ≤0.28 number/mm^2.[69]^11 Signature-based score Immune-associated score Exhausted T cells and cytotoxic T lymphocytes (CTL) were calculated using the single-sample gene set enrichment algorithm from the IOBR package. The tumor immune dysfunction and exclusion (TIDE) score was applied to assess the likelihood of tumor immune evasion in the bulk RNA sequencing profile of tumor samples and predicted the response to immune checkpoint blockade ([70]http://tide.dfci.harvard.edu/). TLS-associated score The TLS-50, TLS-9, and 12-chemokine signatures were used to evaluate TLS abundance based on the mean gene expression value.[71]^21 The CXCL9/10 score, identified as the mean expression level of CXCL9 and CXCL10, was used to evaluate the level of pTLS infiltration. Spatial transcriptome analysis Data processing and quality control Four patients with primary HCC with leading-edge sections ([72]http://lifeome.net/supp/livercancer-st/data.htm)[73]^21 were enrolled in the subsequent spatial transcriptome analysis: (1) the H&E slides of HCC encompassed both tumor and peritumoral areas and (2) the peritumoral area exhibited TLS infiltration. Briefly, spots isolated from the main tissue sections were isolated, and genes expressed in fewer than three spots and mitochondrial/ribosomal-associated genes were filtered.[74]^21 Integration and dimensionality reduction analysis Following quality control, the SCTransform function was applied to normalize the spatial transcriptome and standard data, and identify variable gene features. Subsequently, principal component analysis was performed using the RunPCA function based on 3,000 variable features.[75]^21 To reduce batch effects among patients, we iteratively corrected the spots’ low-dimensional principal components representation using the RunHarmony function, and unsupervised shared-nearest-neighbor-based Uniform Manifold Approximation and Projection clusters were applied to the corrected principal component matrix. Signature-based cell type scoring On average, 14 cells were assessed at each spot.[76]^21 We calculated the log-transformed normalization expression value for each cell type according to both the cell type[77]^21 and xCell[78]^22 signatures. In addition, differentially expressed genes (DEGs) were identified between different cell-type groups and pTLS density groups using the FindAllMarkers function. The non-parametric Wilcoxon rank-sum test and false discovery rate were applied to assess statistical significance and control for multiple comparisons. Pathway enrichment analysis The adjusted p value threshold of <0.05 was used to identify DEGs. To investigate the potential biological function of DEGs, Gene Set Enrichment Analysis (GSEA) was conducted based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome gene sets using the ClusterProfiler R package. Metabolic activity was measured using the scMetabolism R package. Trajectory analysis and identify key regulators To further explore the developmental trajectories of the pTLS density group, pseudotime analysis was applied on specific subsets using the Monocle package (V.2.10.1). Briefly, the fold change of DEGs>0.2 were used in the trajectory analysis after applying the DDRTree and OrderCell functions. Subsequently, the potential pedigree differentiation trajectories were inferred using the “differentialGeneTest” function. Key regulators associated with high pTLS density were identified through intersection analysis of significant DEGs related to pTLS density development and those upregulated in the immune cells infiltration region (IIR). Immunohistochemistry analysis HCC tissues randomly obtained from Xiangya Hospital were fixed in formaldehyde and embedded in paraffin. Owing to incomplete tissue sections from 12 patients, 88 patients with HCC derived from the training set were enrolled in the pathological staining analysis. Continuous 5 mm thick tissue sections were used for the relative quantification analysis of CD3 (1:200 dilution, Huabio, HA720082), CD20 (1:1,000 dilution, Proteintech, 60271–1-lg), CXCL9 (1:200 dilution, Proteintech, 22355–1-AP), and CXCL10 (1:500 dilution, Proteintech, 10937–1-AP) in the pTLS area. Briefly, the process involved dewaxed the tissues in xylene, antigen retrieval in 10 mM citrate buffer, blocking with goat serum, and incubation with the primary antibody at 4°C overnight, followed by incubation with the secondary antibodies (1:200 dilution, Servicebio, GB23301 and GB23303) at 37°C for 45 min. Diaminobenzidine staining and hematoxylin counterstaining of the nuclei were performed. Relative quantification of immunohistochemistry (IHC) targets was then performed using the H-score, which correlates with both the proportion of positively stained cells and the intensity of staining for specific proteins. MR images acquisition MR images were acquired from two scanner manufacturers (1.5/3.0 Tesla GE Health and Siemens Health) with various acquisition settings ([79]online supplemental table 1). Considering the accuracy, stability and generalization of the radiomics models, we selected the specific sequences (T2-weighted imaging (T2WI), T1-weighted arterial phase (AP), and portal vein phase (PVP)) for our study, consistent with previous studies.[80]23,[81]25 Before manual segmentation, all MR images underwent preprocessing (N4 filed bias correction, image standardizations, and voxel resampling) ([82]online supplemental appendix 1). Radiomic analysis Radiomic analysis was conducted according to the Radiomics Quality Score reporting guidelines[83]^26 ([84]online supplemental table 2). Tumor segmentation and feature extraction MR images of each patient were reviewed by two radiologists with 6 and 15 years of experience in abdominal imaging, respectively. Any discrepancies in lesion boundaries were resolved through consensus-based discussion. The tumor margin was manually delineated and the intratumoral area was calculated using ITK-snap software (V.3.8.0, [85]http://www.itksnap.org/). A 5 mm (depending on pixel size) peritumor region was generated by automatically expanding the lesion boundaries, and the intratumoral area was removed. Intratumor and peritumor feature extraction were performed using PyRadiomics software (V.3.0.1, [86]https://github.com/AIM-Harvard/pyradiomics). To ensure reliability, 40 patients were randomly selected, to and the same segmentation was repeated within 1 month to calculate the intraclass and interclass correlation coefficients (ICC) as described previously.[87]^26 Details are provided in [88]online supplemental appendix 2-4. Features selection Variable standardization was performed across the above data sets. To select stable features, the original radiomic features extracted from intratumoral and peritumoral regions were reduced to 2024 (T2WI: 951, AP: 347, and PVP: 726) and 1903 (T2WI: 746, AP: 582, and PVP: 575) based on intra-ICC and inter-ICC values >0.75 ([89]online supplemental appendix 4). Minimum redundancy and maximum relevance were used to further screen the 60 features. Finally, the least absolute shrinkage and selection operator was applied to extract 29 non-zero weighted peritumoral and 30 intratumoral features based on N-fold (Intra: eightfold, Peri: 10-fold) cross-validation ([90]online supplemental figures 1,2) and bootstrap resampling was performed (1,000 times) for radiomic model construction ([91]online supplemental figure 3). Radiomic model construction Three models were built using selected radiomic features ([92]online supplemental table 3): (1) the Intra model, which included only intratumoral features; (2) the Peri model, which included only peritumoral features; and (3) the merged (mPI) model, which merged peritumoral and intratumor radiomic features. The model formula and details are presented in [93]online supplemental appendix 5. Model evaluation The bootstrapping algorithm (1,000 bootstrap intervals) was used to calculate a 95% CI for measuring the model’s performance. Based on the model predictive radiomic-score, patients were classified as having a high or low pTLS density using the maximum Youden Index as the optimal cut-off value. The predictive performance of each model was assessed based on area under the curve (AUC), sensitivity, specificity, accuracy, and F[1] score. A calibration curve was used to assess the goodness-of-fit. Decision curve analysis (DCA) was performed to evaluate the clinical utility. To improve model interpretability, GSEA was used on DEGs between the model-based groups in a subset of the TCGA-LIHC cohort (n=14) with MRI data. Statistical analysis Imaging preprocessing and feature extraction were conducted using Python software (V.3.7; Python Software Foundation). All data cleaning, feature selection, model construction, visualization, and statistical analysis were performed using R software (V.3.4.1; R Foundation for Statistical Computing). Quantitative statistics are presented as the mean (SD) or median (IQR). The Wilcoxon rank-sum test or Kruskal-Wallis test was used for continuous variables with appropriate multiple corrections (Benjamini-Hochberg or Holm-Bonferroni). The χ^2 test or Fisher’s exact test was used for categorical variables. Spearman’s correlation analysis was performed to assess the association between the features. The DeLong test was used to compare the predictive performance of the models. Further, the net reclassification index, integrated discrimination improvement, and likelihood ratio test were employed to assess the improvement between the single-modality and optimal models. Kaplan-Meier (K-M) curves were compared using the log-rank test. Multivariate Cox regression analysis was used to calculate the 95% CIs of HRs to ensure the prognostic value of different variables. The significance threshold was set at p<0.05 (two-sided). Results Baseline characteristics of four in-house HCC cohorts To systematically evaluate the clinical applicability of high pTLS density in HCC, we collected data derived from four international multicenter cohorts (n=952) and four Chinese in-house cohorts (n=701) and immunotherapy in solid tumor cohorts (HCC, n=65; NSCLC, n=28; pan-cancer, n=100; and BLCA, n=298). Subsequently, we integrated the spatial transcriptomics data set (high pTLS density: HCC2 and HCC3, low pTLS density: HCC1 and HCC4) to determine the relevance of high pTLS density as a biomarker of HCC and further verified this using IHC staining. MRI data from patients with HCC in the in-house cohorts were used to develop a non-invasive radiomic model for predicting high pTLS density ([94]figure 1A). Figure 1. Impact of peritumoral tertiary lymphoid structure abundances on prognosis and therapeutic response in hepatocellular carcinoma. (A) Schematic overview of the study. (B and C) Kaplan-Meier analysis of overall survival and relapse-free survival between pTLS density group in patients with HCC who underwent surgery enrolled in the TCGA-LIHC and in-house HCC cohorts. (D) The proportion of immunotherapy response subpopulations in the pTLS^high group compared with pTLS^low groups in outcome cohort (n=41). (E) Distribution of immune-related genes and their relative signatures between pTLS density groups in the TCGA-LIHC cohort. P values were calculated using the Wilcoxon test with Benjamini-Hochberg correction. (F) Tumor immune dysfunction and exclusion score comparing the pTLS^high (n=107) with pTLS^low groups (n=256) in the TCGA-LIHC cohort. BLCA, urothelial bladder carcinoma; HCC, hepatocellular carcinoma; mPI, merged peritumoral and intratumoral radiomic model; NSCLC, non-small cell lung cancer; NR, non-response; PD, progress disease; PD-1, programmed cell death protein-1; PD-L1, programmed death-ligand 1; PR, partial response; pTLS^high, high peritumoral tertiary lymphoid structure density; pTLS^low, low peritumoral tertiary lymphoid structure density; pTLS, peritumoral tertiary lymphoid structures; R, response; SD, stable disease; TCGA-LIHC, The Cancer Genome Atlas-Liver HCC. [95]Figure 1 [96]Open in a new tab Ultimately, the in-house cohorts included 701 patients with HCC, comprising 307 patients in the training data set (64 (20.8%) high pTLS density, mean age, 52 years (IQR 45–62); 53 men (82.8%)), 76 patients in the temporal set (21 (27.6%) high pTLS density, mean age, 55 years (IQR 48–64); 41 men (74.5%)), 277 patients in the external validation set (91 (32.8%) high pTLS density, mean age, 54 years (IQR 48–65); 83 men (91.2%)), and 41 patients in the outcome cohort (16 (39.0%) high pTLS density, mean age, 51 years (IQR 43–54); 14 men (87.5%)) ([97]table 1). Table 1. Baseline patient characteristics. Training set (n=307) Temporal set (n=76) External validation set (n=277) Outcome cohort (n=41) Characteristics pTLS^low (n=243) pTLS^high (n=64) pTLS^low (n=55) pTLS^high (n=21) pTLS^low (n=186) pTLS^high (n=91) pTLS^low (n=25) pTLS^high (n=16) Age (year)[98]^* 52.0(45.0–61.0) 52.0(45.0–62.5) 52.0(48.0–62.0) 55.0(48.0–64.0) 57.0(49.2–65.0) 54.0(48.0–65.5) 58.0(50.0–68.0) 50.5(43.0–53.8) Sex (male) 183 (75.3%) 53 (82.8%) 41 (74.5%) 14 (66.7%) 161 (86.6%) 83 (91.2%) 21 (84.0%) 14 (87.5%) Hepatitis virus infection (present) 175 (72.0%) 51 (79.7%) 37 (67.3%) 14 (66.7%) 159 (85.5%) 80 (87.9%) 21 (84.0%) 14 (87.5%) Cirrhosis (present) 195 (80.2%) 52 (81.2%) 47 (85.5%) 19 (90.5%) 157 (84.4%) 74 (81.3%) 24 (96.0%) 12 (75.0%) AST (>35 U/L) 139 (57.2%) 43 (67.2%) 32 (58.2%) 12 (57.1%) 95 (51.1%) 35 (38.5%) 18 (72.0%) 8 (50.0%) ALT (>40 U/L) 104 (42.8%) 35 (54.7%) 19 (34.5%) 9 (42.9%) 103 (55.4%) 51 (56.0%) 14 (56.0%) 7 (43.8%) TBIL (>20.4 g/L) 87 (35.8%) 19 (29.7%) 25 (45.5%) 7 (33.3%) 53 (28.5%) 33 (36.3%) 4 (16.0%) 3 (18.8%) ALB (>40 g/L) 79 (32.5%) 24 (37.5%) 17 (30.9%) 9 (42.9%) 53 (28.5%) 33 (36.3%) 11 (44.0%) 8 (50.0%) AFP (>400 ng/mL) 83 (34.2%) 19 (29.7%) 20 (36.4%) 10 (47.6%) 53 (28.5%) 31 (34.1%) 12 (48.0%) 5 (31.2%) Tumor number (≥2) 50 (20.6%) 8 (12.5%) 10 (18.2%) 7 (33.3%) 50 (26.9%) 24 (26.4%) 5 (20.0%) 1 (6.25%) Maximum diameter (cm)[99]^* 5.15(3.24–8.13) 5.42(3.57–7.31) 5.89(3.17–8.42) 4.84(3.35–6.84) 5.20(3.32–8.00) 4.50(3.21–6.87) 7.01(4.07–10.8) 5.56(3.29–7.07) [100]Open in a new tab Except where indicated, data are presented as the numbers of patients, with percentages in parentheses. ^* Data are medians, with IQRs in brackets AFPalpha-fetoproteinALTalanine aminotransferaseASTaspartate transaminasepTLSperitumoral tertiary lymphoid structureTBILtotal bilirubin Prognosis and therapeutic response of high pTLS density We further investigated the clinical relevance of high pTLS density by examining its influence on survival and immunotherapy response in patients with HCC. High pTLS density group significantly correlated with improved overall survival (OS) and relapse-free survival (RFS) (both log-rank p<0.05)in the TCGA-LIHC (n=363) and in-house (n=660) cohorts ([101]figure 1B–C). In the dual-center outcome cohort (n=41) ([102]figure 1D), pTLS^high subset levels were higher in treatment responders than in non-responders (43.7%, p=0.03). Moreover, TIDE scores were significantly lower in the pTLS^high group than in the pTLS^low group (p=0.043), suggesting the significance of high pTLS density in identifying patients who are likely to benefit from anti-PD(L)-1-based treatment strategies ([103]figure 1E). In the TCGA cohort, the presence of lymphocyte infiltration (CTLs, B lineage), upregulated PD(L)-1 expression, and reduced exclusive T-cell infiltration were observed in the pTLS^high group, collectively supporting a correlation between favorable clinical outcomes ([104]figure 1F). In addition, GSEA revealed immune inflammation pathway activation in the pTLS^high group (p<0.001) ([105]online supplemental figure 4). These findings demonstrated that high pTLS density status impacts the anticancer microenvironment, although its specific biological function remains elusive. General profiles of spatial transcriptomics To comprehensively understand the spatial heterogeneity and biological function of high pTLS density, we reanalyzed the spatially sequential sections collected from the leading-edge section of four HCC tumors (two high pTLS density and two low pTLS density). As each spot contained several cells, we annotated the spots using unbiased and comprehensive strategies, including signature, cluster, and marker genes according to previous reports.[106]^21 Subsequently, the slides of each patient were categorized into five sections based on histological and cluster results, including tumor section, fiber cord sector, peritumoral hepatocytes, peritumoral IIR and other sections ([107]figure 2A and [108]online supplemental figure 5A). Differential expression marker genes elevated in other regions, such as the tumor region (GPC3, DLK1, SCD), hepatocytes (MT1G, CYP2E1, CYP3A4), and fiber cord sector (COL1A1, COL1A2, ACTA2) were identified (all p<0.001, [109]online supplemental figure 5B). TLS annotation was further validated through histological analysis ([110]figure 2B), the spatial location of CD3/CD20, and TLS signatures, and indicated that the IIR region corresponded to the pTLS region (p<0.001, [111]online supplemental figure 5C). Figure 2. Leading-edge section of spatial transcriptomics in hepatocellular carcinoma. (A and B) H&E-stained tissue sections from the spatial transcriptomic cohorts (n=4). Cell types (A, color spots) and peritumoral tertiary lymphoid structures (B, bold numbers) are annotated and marked. (C) Dot plot showing immune cell signature-based scores in the pTLS^high and pTLS^low groups within the immune-infiltration region. The dot size indicates the percentage expression. P values were calculated using the Wilcoxon test with Holm-Bonferroni correction. (D) TLS-associated gene markers in the pTLS^high group compared with pTLS^low group within the immune-infiltration region. The dot color indicates average expression. P values were calculated using the Wilcoxon test with Holm-Bonferroni correction. (E) Volcano plot showing the upregulated differential genes between the pTLS groups in the immune-infiltration region of patients with HCC. P values were calculated using the Wilcoxon test with Benjamini-Hochberg correction. (F) Pseudotime analysis revealing the cell trajectory of spots in the IIR from the start (blue) to the endpoint (red), colored by states and pseudotime. (G) Developmental trajectory of spots between the pTLS^high (red) and pTLS^low (blue) groups. (H) Venn diagram of the top 50 differentially expressed genes that were upregulated in the IIR and during trajectory analysis. P values were calculated using the Wilcoxon test with Benjamini-Hochberg correction. (I) Developmental trajectory of CXCL9 and CXCL10. (J) Bubble plot showing the mean expression level between IIR and non-IIR sections. The dot size indicates the mean value of the genes. p>0.05; ***p<0.001. Avg.exp, average expression; Bm, mature B cell; Bn, Bn, naïve B cell; GC, germinal center; HCC, hepatocellular carcinoma; ICIs, immune checkpoint blockades; IFN-G, interferon gamma; IIR, immune-infiltration region; ns, non-significant; pTLS^high, high peritumoral tertiary lymphoid structure density; pTLS^low, low peritumoral tertiary lymphoid structure density; pTLS, peritumoral TLS; TLS, tertiary lymphoid structures. [112]Figure 2 [113]Open in a new tab Spatial immune characteristics of high pTLS density Immune cell infiltration difference depends on the pTLS status. The pTLS^high group exhibited more B cell, dendritic cell, macrophage lineages and T central memory cells in the IIR (all p<0.001, [114]figure 2C). Further B lineages subpopulation analysis revealed that the pTLS^high group exhibited more plasma cells (JCHAIN^+/IGHA1^+/IGHG4^+), memory B cells (CD27^+/TNFRSF13B^+/SDC^−), and germinal center (CD38^+/MKI67^+/CXCR4^+) rather than naive B cells (MS4A1^+/IL4R^+) in the IIR (all p<0.05, [115]figure 2D). In addition, immune checkpoint inhibitor (ICI) positive T cells (PDCD1/HAVCR2/LAG3/CD274) and inflammatory factors (HLA-DRA/IDO1/CXCL9) were upregulated in the pTLS^high group (all p<0.05), indicating a potential for improved ICI response rates, consistent with our transcriptomic findings. To further investigate the potential mechanism, differential expression analysis was performed to identify DEGs associated with key signaling pathways. In the volcano plot, inflammatory microenvironment-associated genes (ICAM1/CXCL9/CXCL10), humoral immunity markers (IGHA1/IGHG1-4), and toll-like receptor genes (SPP1/CXCL9) were significantly upregulated in the pTLS^high group (all p<0.001,[116]figure 2E). Furthermore, KEGG and Reactome enrichment analyses revealed that high pTLS density is correlated with enhanced antigen processing and presentation, increased inflammatory immune cell infiltration, and persistent antitumor immune memory establishment (all p<0.05, [117]online supplemental figure 6A). At the metabolic level, the pTLS^low group was enriched in tumorigenesis-associated metabolic pathways, such as glutathione and glycolysis/gluconeogenesis (both p<0.01, [118]online supplemental figure 6B). Collectively, these results suggested that high pTLS density was associated with antitumor-associated inflammatory immune cell infiltration, thereby contributing to the defense against HCC and benefits from ICI-based treatment. CXCL9/10 score was associated with the occurrence of high pTLS density To further identify the key regulator of high pTLS density formation, pseudotime analysis was performed and revealed that pTLS^high status was primarily concentrated at the early stages and the endpoint of branch one in the developmental trajectory of the pTLS^low group ([119]figure 2F–G). To identify DEGs at the early stage, intersection analysis was conducted to screen the top 14 genes between the pseudotime analysis and IIR-related DEGs (top 50, p<0.001, [120]figure 2H). In addition to the TLS-associated genes, CXCL9 and CXCL10 co-accumulation at the starting point of the developmental trajectory was observed, which decreased along the pseudotime trajectory, suggesting that CXCL9 and CXCL10 co-expression may be related to the occurrence of high pTLS density ([121]figure 2I and [122]online supplemental figure 6C,D). Both CXCL9 and CXCL10 were predominantly secreted by dendritic cells and macrophages in areas enriched with immune checkpoint genes[123]^27 28 and co-recruited via LAG3^+/CXCL13^+ T cells in the melanoma, promoting both naïve T-cell infiltration and the inflammatory tumor microenvironment.[124]^29 While various studies have demonstrated that increased CXCL9 or CXCL10 expression correlates with therapeutic response in patients with advanced HCC who received atezolizumab plus bevacizumab,[125]^30 31 few studies have reported their co-accumulation in HCC. Hence, we hypothesized that co-upregulation of CXCL9 and CXCL10 may be the leading cause of the development of high pTLS density in IIR. Compared with when expressed alone, the combined CXCL9/10 score was significantly higher in the IIR in the high pTLS density group (all p<0.001, [126]figure 2J and [127]online supplemental figure 6E), and positively correlated with TLS-associated signatures, 12-chemokine markers, and plasma cells in the IIR ([128]online supplemental figure 6F). These findings suggest that CXCL9 and CXCL10 co-accumulation may promote immune cell infiltration and a pro-inflammatory pTLS region. Predictive role of CXCL9 and CXCL10 in favorable HCC outcomes To further verify the spatial distribution between pTLS density and CXCL9/10 score, patients with HCC from the Xiangya cohort with leading-edge were enrolled. Consistent with the spatial analysis, patients with HCC exhibited in the pTLS^high group have increased infiltration of CXCL9 and CXCL10 in the IIR ([129]figure 3A–B). Furthermore, the H-score analysis indicated that the high pTLS density group exhibited significantly upregulated CXCL9/CXCL10 expression compared with low pTLS density (n=88, p<0.001, [130]figure 3C). Similarly, the TCGA-LIHC cohort exhibited the same results (n=363, p=0.025, [131]figure 3D). Additionally, CXCL9/10 score positively correlated with TLS-50 signature across the four international multicenter HCC cohorts (n=952, Spearman’s coefficient: 0.49–0.55, p<0.001, [132]figure 3E), which was consistent with our spatial analysis. Figure 3. High CXCL9/10 score predicts a favorable prognosis and immune-checkpoint blockades response. (A and B) Representative immunochemical staining of CXCL9, CXCL10, CD3, and CD20. (C and D) Differential distribution of CXCL9/CXCL10 score/H-score between the pTLS density groups in the Xiangya and TCGA-LIHC cohorts. H-score=immunohistochemistry-score (E) correlation between CXCL9/CXCL10 score and TLS-50 signature across the four HCC cohorts. (F and G) Kaplan-Meier survival analysis shows the correlation between CXCL9/10 expression patterns and pTLS density groups in the TCGA-LIHC and Xiangya cohorts. (H) Kaplan-Meier survival analysis showings the survival difference between CXCL9/10 score-based optimal cut-off values in [133]GSE76427, TCGA-LIHC, and [134]GSE14520 cohorts. (I) Boxplot showing the CXCL9/10 score expression levels between non-responders and responders in NSCLC (n=27), HCC (n=21), pan-cancer (n=100), and BLCA (n=298) cohort. (J) Prognostic differences between the CXCL9/10 score-derived group across three immunotherapy cohorts. (K) Receiver characteristic curve analysis showing the predictive applicability of CXCL9/10 score for immunotherapeutic response across four cohorts. BLCA, bladder cancer; CR, complete response; HCC, hepatocellular carcinoma; ICHC-LIRI, International Cancer Genome Consortium-Liver Cancer RIKEN; NSCLC, non-small cell lung cancer; PD(L)1, programmed cell death protein or receptor 1; PR, partial response; pTLS, peritumoral TLS; SD, stable disease; TCGA-LIHC, The Cancer Genome Atlas-Liver HCC; TLS, tertiary lymphoid structures. [135]Figure 3 [136]Open in a new tab Subsequently, we evaluated the prognostic and immunotherapeutic effects of CXCL9/10 score. K-M survival analysis revealed that the CXCL9/10^highpTLS^high group had significantly prolonged mST for RFS (or disease-free survival (DFS)) (p<0.05) and OS (p<0.05) compared with the other subtypes in both the TCGA-LIHC and Xiangya cohorts ([137]figure 3F–G). Compared with a single gene, a high CXCL9/10 score was associated with favorable mST for OS or RFS (or DFS) across the four HCC cohorts (all p<0.05, [138]figure 3H and [139]online supplemental figure 7). Further treatment evaluation demonstrated that anti-PD(L)-1 treatment responders had a higher CXCL9/10 score in various solid tumors cohorts (n=446, all p<0.05, [140]figure 3I) than that of non-responders. The survival difference based on the CXCL9/10 score was also verified across three immunotherapeutic cohorts with prognostic information (all p<0.05, [141]figure 3J). Among them, the CXCL9/10 score demonstrated the best predictive ability in NSCLC, followed by HCC, and the worst in BLCA (AUC: 0.61–0.85, [142]figure 3K). These results suggest that CXCL9 and CXCL10 may act synergistically to promote high pTLS density formation, thereby affecting prognosis and anti-PD(L)-1-based therapy in patients with HCC by promoting the pro-inflammatory tumor microenvironment. Radiomic model predicts high pTLS density in HCC For the non-invasive detection of the pTLS status of patients with HCC, we developed a classifier based on quantitative MR radiomic features ([143]figure 4A). First, we used patients with HCC who underwent surgery (n=660) in the in-house cohorts and divided them into a training (n=307), temporal test (n=76), and external validation (n=277) set for model construction and validation ([144]figure 4B). Subsequently, 30 intratumor and 29 peritumor radiomic features were used to establish three logistic linear models (Intra, Peri, and mPI) to predict pTLS abundance. Compared with the individual single-modality models, the mPI model showed a significant improvement in predictive performance according to the likelihood ratio analysis ([145]online supplemental table 4, all p<0.001), Integrated Discrimination Improvement Index ([146]online supplemental table 5, all p<0.001), and Net Reclassification Index ([147]online supplemental table 6, all p<0.05). Regarding the discrimination ability and generalizability, the mPI model exhibited superior performance for identifying high pTLS density, with an AUC of 0.91 (95% CI (0.87, 0.94)) than those of the individual single-modality models in the external validation set (DeLong test, both p<0.05) ([148]table 2 and [149]online supplemental table 4). Additionally, for model stability under different conditions, no significant difference was observed in the AUC between the following intergroup comparisons: training versus external set (DeLong test, p=0.61), external versus temporal set (DeLong test, p=0.89), and training versus temporal set (DeLong test, p=0.87) ([150]online supplemental table 7). Similarly, the DCA indicated that the mPI model had a greater net benefit for detecting high pTLS density than those of the individual single-modality models ([151]online supplemental figure 8A-C). For goodness-of-fit ability, the calibration curve showed stable performance in three radiomic models ([152]online supplemental table 8, [153]online supplemental figure 8D-F). For biological interpretability, the predicted pTLS^high group exhibited significant TLS-associated pathway activation in the TCGA-LIHC cohorts with enhanced-contrast MRI data (n=14, false discovery rate <0.05) ([154]online supplemental table 9). Figure 4. Radiomic signature predicts high pTLS density. (A) Workflow of the radiomic-based model. (B) Flowchart of patient inclusion and exclusion among the multicenter data sets. (C–E) Area under the receiver operating characteristic curve (AUC) analysis for predicting high pTLS density between the mPI, Peri, and Intra models in the training (C) and temporal (D) and external validation (E) sets. (F–H) Survival differences for relapse-free survival (upper) and overall survival (lower) between the predicted high and low pTLS density groups using the mPI model across three data sets. (I) AUC analysis for predicting therapy response based on mPI model scores in the outcome cohort. (J and K) Distribution of standard mPI scores between the PR, n=10 and PD/SD, n=31 groups. P values were calculated using the Wilcoxon test. AP, arterial phase of T1-weighed image; Intra, intratumoral; mPI, merged Peritumoral and Intratumoral MRI rad-score; Peri, peritumoral; PD, progress disease; PR, partial response; PVP, portal vein phase of T1-weighted image; pTLS^high and pTLS^low, high and low pTLS density; pTLS, peritumoral tertiary lymphoid structures; ROC, Receiver operating characteristic curve; SD, table disease; Radscore, radiomic-derived model score; T2WI, T2 weighted images. [155]Figure 4 [156]Open in a new tab Table 2. Predictive performance of radiomic model for predicting high pTLS density. Training set Temporal set External validation set Parameter (n=307) (n=76) (n=227) mPI model  AUC 0.92 (0.89, 0.96)[157]^* 0.91 (0.83, 0.99)[158]^* 0.91 (0.87, 0.94)[159]^*  Sensitivity (%) 88 (77, 97) (56/64) 86 (67, 100) (18/21) 91 (86, 97) (83/91)  Specificity (%) 84 (75, 92) (205/243) 95 (84, 100) (52/55) 84 (78, 89) (156/186)  Accuracy (%) 85 (261/307) 91 (70/76) 86 (239/277)  F[1]-score 0.71 0.84 0.81 Intra model  AUC 0.87 (0.82, 0.91) 0.82 (0.71, 0.93) 0.84 (0.79, 0.89)  Sensitivity (%) 84 (66, 97) (54/64) 90 (48, 100) (19/21) 70 (58, 91) (64/91)  Specificity (%) 77 (61, 92) (188/243) 67 (49, 100) (37/55) 86 (59, 94) (160/186)  Accuracy (%) 79 (242/307) 75 (56/76) 80 (224/277)  F[1]-score 0.61 0.66 0.69 Peri model  AUC 0.85 (0.80, 0.90) 0.77 (0.64, 0.90) 0.76 (0.70, 0.82)  Sensitivity (%) 81 (69, 92) (52/64) 81 (38, 100) (17/21) 64 (48, 88) (58/91)  Specificity (%) 78 (67, 89) (190/243) 71 (49, 100) (39/55) 81 (54, 91) (151/186)  Accuracy (%) 79 (242/307) 75 (56/76) 74 (209/277)  F[1]-score 0.61 0.64 0.63 [160]Open in a new tab AUC, sensitivity, specificity, and accuracy in the three data sets were calculated with cut-off values that maximized the Youden Index in the training data set. The specific cut-off values were calculated using the Youden Index based on bootstrapping (repeated 1000 times): −1.754 for the mPI model, 0.145 for the Intra model and 0.208 for the Peri model. Except where indicated, numbers in parentheses are presented as 95% CIs and numbers in brackets present the numbers of patients. ^* Statistical significance calculated using the DeLong test of AUC between single-modality models and the mPI model. ACUarea under curveIntraintratumoralmPI modelmerged Peri and Intra MRI radiomic modelPeriperitumoral Radiomic model predicts prognosis and immunotherapy response in HCC Patients with HCC who underwent surgery were stratified into the predicted pTLS^high and pTLS^low groups based on the optimal cut-off value using the training set to evaluate potential clinical utility. The predicted pTLS^high group exhibited prolonged mST for OS (HR 0.20–0.44, p=0.005 to <0.001) and RFS (HR 0.25–0.73, p=0.015 to <0.001) compared with predicted pTLS^low group across three data sets ([161]figure 4F–H). Regarding other prognostic factors (eg, resection margin and vascular invasion status[162]^32), multivariable Cox regression analysis revealed that the mPI model was an independent predictor using RFS (HR=0.56, (95% CI 0.46, 0.69), p<0.001) and OS (HR=0.51, (95% CI 0.40, 0.66), p<0.001) ([163]online supplemental table 10). Additionally, the radiomics-based score was positively correlated with the CXCL9/10 H-score in the Xiangya cohort ([164]online supplemental figure 9, p<0.001, n=88). Next, we investigated the predictive role of radiomics-based scores in predicting immunotherapy response in HCC. Treatment response assessment revealed that among the 41 included patients with HCC who received immunotherapy plus anti-angiogenic therapy, 8 (19%) responders exhibited a partial response and 33 (81%) non-responders demonstrated stable or progressive disease ([165]online supplemental table 11, [166]online supplemental figure 10). Similarly, patients with advanced HCC exhibiting high pTLS density experienced a longer mST for OS and progression-free survival ([167]online supplemental figure 11, both p<0.05). The mPI score was found to be a good predictor of treatment response ([168]figure 4I, AUC=0.835). The responders exhibited increased baseline mPI model scores compared with those of non-responders (p=0.009) ([169]figure 4J–K), while no significant difference was observed in the Intra and Peri model scores ([170]online supplemental figure 12). Overall, the optimal radiomics-based model score was identified as an independent predictor of prognosis, thereby underscoring the predictive value of the ICI-based anti-angiogenic therapy response in HCC. Discussion Previous studies have confirmed that the presence of TLS is associated with favorable prognosis and enhanced immunotherapy-based response across various solid tumors.[171]^4 6 18 33 In this study, we validated the clinical utility of high pTLS density across multicenter cohorts, identified a potential biomarker using multiomics resolution and proposed an enhanced-contrast MRI-based radiomic model to non-invasively predict high pTLS density in HCC. Consistent with previous study,[172]^11 we confirmed the prognostic value of pTLS classification in these cohorts. Notably, in patients with advanced HCC who received anti-PD(L)-1 and antiangiogenic treatment in dual-center cohorts, those with high pTLS density showed better responses. Similarly, another study found that increased TLS density was positively correlated with major or complete pathological response in 19 patients with HCC who underwent neoadjuvant ICI.[173]^10 Matched pathway analysis reveals consistent results with our study regarding T and B-cell activation, suggesting a relationship between preferential immunotherapy response and increased pTLS density. Furthermore, we observed that high pTLS density correlated with B-cell profile and antitumor signaling pathway, which aligned with previous reports.[174]^4 7 15 34 Interestingly, unsaturated fatty acid and glycolytic metabolism were associated with low pTLS density, suggesting that tumorigenesis-associated metabolic processes may affect TLS formation. Therefore, our study provides evidence of an association between pTLS density and anti-PD(L)-1-based combination treatment response. Solid tumors and TLS have been associated with chemokine infiltration, including CCL19, CCL21, and CXCL13.[175]^4 7 15 34 In the present study, we further identified potential chemokines (CXCL9 and CXCL10) that promote high pTLS density through spatial analysis and multicenter verification. Elevated levels of CXCL9 and CXCL10 have been confirmed as predictors of immunotherapy response in HCC.[176]^27 30 Furthermore, compared with other chemokines, CXCL9 and CXCL10 were significantly enriched in exhausted CD8^+ T-cell milieus in melanoma,[177]^35 suggesting substantial local antitumor activity. Notably, compared with individual expression, CXCL9 and CXCL10 co-expression was consistent with the differential infiltration of exhausted markers, LAG3 and PD(L)-1, prognosis, and the efficacy of immune therapy between high/low pTLS density groups, emphasizing the synergistic effect of these two cytokines. However, further research is needed to verify the specific mechanism of CXCL9 and CXCL10 co-expression in maintaining TLS function and promoting an effective antitumor immune response in patients with HCC and high pTLS density. Radiological-based models have emerged as a non-invasive method to predict the presence of TLS in various solid tumors.[178]78 13,[179]15 36 Compared with a previous CT-based model (mean AUC of 0.75) from a single-center study,[180]^8 our MRI-based radiomic classifier demonstrated superior performance (AUC of 0.91, 95% CI (0.87, 0.94)) in identifying TLS abundance across multicenter HCC data sets. Moreover, we found that the combined intratumoral and peritumoral radiomic feature model performed better than both of the single-modality models for predicting pTLS abundance in the external validation set (AUC range 0.76–0.84, both p<0.001). In another MRI-based study,[181]^14 a multimodality MRI-based radiomic model demonstrated significant model improvement compared with single-modality models, supporting the added value of multimodality radiomic classifiers. The biological interpretability of radiomics is considered a scientific challenge.[182]^26 37 To date, previous radiomic models for predicting TLS in solid tumors have lacked biological interpretability.[183]813,[184]15 36 Recently, Li et al’s study has elucidated the potential biological mechanisms underlying a predictive CT-based radiomic model of intratumoral TLS based on single-cell results, revealing a correlation with mammalian target of rapamycinm (mTOR) signaling pathway activation in TLS.[185]^7 Unlike that study, we identified potential combined drivers of pTLS density formation through spatial dimensions, namely the CXCL9/10 score. Our MRI-based radiomic model was significantly correlated with the CXCL9/10 score and the activation of specific TLS-associated pathways (eg, mTOR signaling, B and T-cell activation, and major histocompatibility complex II antigen presentation). These findings provide a solid biological basis for using MRI-based radiomic classifiers to predict pTLS, facilitating future translational applications. Additionally, our study found that the radiomic model for predicting high pTLS density was associated with the objective response rate of anti-PD(L)-1-based combination therapy. However, translation into clinical practice remains challenging, particularly with regard to determining the optimal cut-off values. Further prospective studies involving a broader spectrum of patients with HCC are needed for validation. Our findings are promising; however, there are some limitations. First, as this data set was a retrospectively recruited data set, potential selection bias may exist. Despite efforts to mimic a prospective cohort using a temporal sampling strategy, verification in a prospective trial is warranted. Second, the manual delineation of the area of interest in this study was both time-consuming and labor-intensive. Nonetheless, manual segmentation remains the standard owing to its high precision in MRI radiomics for HCC research. Future advancements should focus on developing an automated segmentation system based on image registration across multiple phases. Third, most patients were diagnosed with virus-related HCC. Therefore, additional validation is essential before extending the applicability of the model to individuals with non-alcoholic steatohepatitis-associated HCC. Fourth, the use of unmatched multiomics data limited the depth and breadth of the research. Finally, various anti-PD(L)-1 antibodies were employed within the study cohort, potentially introducing internal variability given the constrained sample size. However, the combination of angiogenesis inhibitors with distinct immunotherapy antibodies reflects indicative of in-house clinical practice. Conclusion This study evaluated the clinical value of high pTLS density in terms of prognosis and immunotherapeutic response in patients with HCC, and identified spatial expression patterns and potential regulators of pTLS density. Moreover, we constructed a radiomic classifier for the non-invasive prediction of high pTLS density and found a correlation with favorable prognosis and ICI-based combined therapy response. Our novel approach has the potential to optimize precision medicine, through further prospective evidence is needed. supplementary material online supplemental table 1 [186]jitc-12-12-s001.docx^ (52.1KB, docx) DOI: 10.1136/jitc-2024-009879 online supplemental figure 1 [187]jitc-12-12-s002.docx^ (3.3MB, docx) DOI: 10.1136/jitc-2024-009879 online supplemental file 1 [188]jitc-12-12-s003.docx^ (29.5KB, docx) DOI: 10.1136/jitc-2024-009879 Acknowledgements