Noninvasive evaluation of neutrophil extracellular traps signature predicts clinical outcomes and immunotherapy response in hepatocellular carcinoma

Background Neutrophil extracellular traps (NETs) have been shown to play a pivotal role in promoting metastasis and immune escape in hepatocellular carcinoma (HCC). Therefore, noninvasive tests to detect the formation of NETs in tumors can have significant implications for the treatment and prognoses of patients. Here, we sought to develop and validate a computed tomography (CT)-based radiomics model to predict the gene expression profiles that regulate the formation of NETs in HCC. Methods This study included 1133 HCC patients from five retrospective cohorts. Based on the mRNA expression levels of 69 biomarkers correlated with NET formation, a 6-gene score (NETs score, NETS) was constructed in cohort 1 from TCIA database (n=52) and validated in cohort 2 (n=232) from ICGC database and cohort 3 (n=365) from TCGA database. And then based on the radiomics features of CT images, a radiomics signature (RNETS) was developed in cohort 1 to predict NETS status (high- or low-NETS). We further employed two cohorts from Nanfang Hospital (Guangzhou, China) to evaluate the predictive power of RNETS in predicting prognosis in cohort 4 (n=347) and the responses to PD-1 inhibitor of HCC patients in cohort 5 (n=137). Results For NETS, in cohort 1, the area under the curve (AUC) values predicting 1, 2, and 3-year overall survival (OS) were 0.836, 0.879, and 0.902, respectively. The low-NETS was associated with better survival and higher levels of immune cell infiltration. The RNETS yielded an AUC value of 0.853 in distinguishing between high-NETS or low-NETS and patients with low-RNETS were associated with significantly longer survival time in cohort 1 (P<0.001). Notably, the RNETS was competent in predicting disease-free survival (DFS) and OS in cohort 4 (P<0.001). In cohort 5, the RNETS was found to be an independent risk factor for progression-free survival (PFS) (P<0.001). In addition, the objective response rate of HCC patients treated with PD-1 inhibitor was significantly higher in the low-RNETS group (27.8%) than in the high-RNETS group (10.8%). Conclusions This study revealed that RNETS as a radiomics biomarker could effectively predict prognosis and immunotherapy response in HCC patients.


Introduction
Hepatocellular carcinoma (HCC), the most commonly occurring primary liver cancer, remains a global health challenge and is the third leading cause of cancer-related death worldwide (1). The tumor immune microenvironment (TIME), which is composed of various tumor-infiltrating immune cells, is thought to greatly affect tumor progression and response to immunotherapy (2). Currently, blocking the signaling of the programmed cell death receptor-1 (PD-1) and programmed cell death receptor ligand-1 (PD-L1) pathways with monoclonal antibodies is a new immunotherapeutic approach for treating advanced HCC cases. However, the response rates to this treatment have been low and the overall prognosis for HCC patients does not seem to improve with this treatment, probably because of the TIME heterogeneity (3)(4)(5). Therefore, tools that can accurately determine the TIME status could be invaluable in guiding treatment for HCC patients.
Neutrophils that infiltrate tumors, also called tumor-associated neutrophils (TANs), often play pro-tumorigenic roles and support tumor progression in response to microenvironmental cues released by tumor and stromal cells (6,7). Neutrophil extracellular traps (NETs), which are extracellular web-like structures of DNA chromatin complexes extruded from dying neutrophils, were once thought to mainly function as snares that caught and killed harmful microorganisms (8). There is increasing evidence to prove that in cancers, NETs play important role in promoting metastasis (9)(10)(11). Recent analyses suggested that studies on NETs are becoming increasingly momentous in research on liver cancer. For example, one study indicates that the elimination of NETs may slow the progression of steatohepatitis-related liver cancer (12). In addition, the formation of NETs is known to trigger pro-tumorigenic inflammatory responses and fuel HCC metastasis (13). Furthermore, NETs have been reported to protect cancer cells from being attacked by the immune system and reduce the effectiveness of immunotherapy (14)(15)(16). Thus, scoring model that can evaluate the formation potential of NETs in tumors can be useful as biomarkers for predicting survival and in estimating if immunotherapy could be of benefit to HCC patients.
By using the expression profiles of genes associated with NETs formation in biopsy specimens, we are able to directly measure the abundance of NETs in tumor tissues. Recent studies have shown that the integration between machine learning algorithms and NETs-related gene signatures own the potential to be reliable biomarkers in predicting prognosis and the response to immunotherapy in several types of malignant tumors (17,18). However, this method has some major limitations, including the requirement for an invasive procedure and a forbiddingly high cost for patients. The emergence of radiomics, which uses quantitative medical imaging features, maybe a viable alternative to the method mentioned previously (19,20). Computed tomography (CT) is now in standard use for diagnosing, staging and monitoring therapeutic efficacy in liver tumors. CT-based radiomics have shown that it is capable of achieving successful assessment and predictive ability in prognostic aspect for HCC patients (21)(22)(23). Given the increasing use of immunotherapy in advanced HCC cases, the role of CT radiomics in characterizing TIME necessitates more exploration and amelioration. Some research has shown that certain radiological features of tumors are closely related to specific gene expression profiles; this provides a unique opportunity for developing non-invasive complementary approaches to detect NETs in the TIME (24,25).
This study has three objectives. One, to identify a NETs-related gene signature (NETs score or NETS) that correlated with patients' prognosis and could differentiate 'cold' or 'hot' TIME. Two, to establish a link between the NETs score and tumor radiomics features obtained from contrast-enhanced CT. And three, to build a novel radiomics model (radiomics NETs score or RNETS) to predict survival and to identify HCC patients for whom immunotherapy would be beneficial.

Study design and data collection
This study retrospectively included five cohorts and the flowchart in Figure 1 describes the entire study design. The HCC patients with RNA sequence data, baseline CT images and the corresponding clinical data from The Cancer Imaging Archive (TCIA, n=52) were collected as training cohort (Cohort 1) to firstly develop a prognostic NETs-related gene signature (NETs score, NETS), and then construct a CT-based radiomics model (radiomics NETs score, RNETS) that correlated with expression patterns of NETS by using machine learning algorithms. Two external validation cohorts with RNA sequence data and complete follow-up data were downloaded from International Cancer Genome Consortium (ICGC, n=232, Cohort 2) and The Cancer Genome Atlas (TCGA, n=365, Cohort 3) to validate the prognostic value of NETS. Besides, HCC patients with baseline CT information from Nanfang Hospital (Guangzhou, China) who underwent curative resection between January 2013 and October 2018 (Cohort 4, n=347) were used to evaluate the predictive value of RNETS in predicting the risk of postoperative recurrence, while those receiving immunotherapy with unresectable HCC between January 2019 and October 2021 (cohort 5, n=137) from Nanfang Hospital were used to evaluate the predictive power of RNETS in predicting the response to immunotherapy.
For patients without previous treatment, CT imaging data at baseline were collected for RNETS establishment and evaluation. Besides, for patients who were refractory to the standard first-line therapy (loco-regional treatments or targeted therapies) and received immunotherapy as second-line treatment, imaging data for analysis were collected based on the inclusion criteria as follows: 1. available CT imaging data within 30 days prior to PD-1 inhibitors therapy; 2. with ≥ 1 measurable target lesion. Follow-up data were collected from the initiation of treatment until 31 October 2022. Disease-free survival (DFS) was defined as the time from surgery to tumor recurrence at any site or death from any cause. Progressionfree survival (PFS) was set to the time period from the execution of immunotherapy to disease progression or death from any cause. The responses to PD-1 inhibitors were defined as partial response (PR), complete response (CR), stable disease (SD), or progressive disease (PD), and evaluated by CT or magnetic resonance imaging (MRI) at 3 months and 6 months according to the RECIST version 1.1 (26).

Establishment and evaluation of NETs score
A total of 69 NETs-related genes were collected from previous studies (8,27) and are shown in Table S1. Cox proportion hazards regression with the least absolute shrinkage and selection operator (LASSO) was applied to tune the optimal value of the penalty parameter (l) and derive the genes with the highest predictive value from training cohort with non-zero coefficients identified. The calculation formula for NETs score is presented as follows: Where n, P, and b denoted the number of selected genes, the normalized gene expression level, and the coefficient index derived from LASSO regression, respectively. The expression level was obtained by log2-transformed FPKM+1 of each gene. The 'survminer' package was employed to pick out the cut-off point for the NETs score, according to which, patients were divided into two groups (high-NETS and low-NETS). Subsequently, Kaplan-Meier curves were plotted and the log-rank method was employed to compare the survival distributions of the two groups. Utilizing 'pROC' R package and 'timeROC' R package to output receiver operating characteristic (ROC) curves and evaluate the predictive value according to the area under the curves (AUC).

Immune infiltration assessment
The CIBERSORT deconvolution algorithm was used to quantify the infiltration levels of 22 classes of immune cells in the high-and low-NETS groups (28). In this step, we analyzed the data of normalized gene expression values combined with the LM22 signature and performed 1,000 permutations in R studio (v4.1.0). ESTIMATE was used to calculate different scores to describe the immune Flowchart of overall study.
Xin et al. 10.3389/fimmu.2023.1134521 Frontiers in Immunology frontiersin.org microenvironment, such as the infiltration level of immune cells (ImmuneScore), and stromal content (StromalScore) (29). Box plots were generated based on the scores obtained by the two NETS groups.

Functional pathway enrichment analysis
Differentially expressed genes (DEGs) between the high-and low-NETS groups were obtained by using the empirical Bayesian approach via the 'limma' package. The significance criteria for selecting the DEGs were set as the false discovery rate (FDR)< 0.05 and the absolute value of the log2 (fold change) was set to ≥ 1.5. We then performed Gene Set Enrichment Analysis (GSEA) using the 'GSEAbase' R package to recognize the hallmark pathways in low-and high-NETS groups respectively.

CT image segmentation and radiomics feature extraction
The CT images from the portal venous phase were fetched for extracting image features. We made use of the ITK-SNAP software (v4.0) (www.itksnap.org) to perform the delineation for the regions of interest (ROIs), in which two radiologists (Reader 1 and Reader 2, both with > 5 years of experience) manually segment the entire tumor on each axial slice. The CT images were resampled to a voxel size of 1×1×1 mm 3 to standardize the voxel spacing. Voxel intensity values were discretized by a fixed bin width of 25 hounsfield unit (HU) to reduce image noise and normalize intensities; wavelet filtering was employed to all the CT series.
Subsequently, 1316 radiomics features were extracted from each ROI using the Pyradiomics Python package (version 3.0). Seven types of features were included: 1) Shape; 2) First Order Statistics; 3) Gray Level Co-occurrence Matrix; 4) Gray Level Run Length Matrix; 5) Gray Level Size Zone Matrix; 6) Neighboring Gray Tone Difference Matrix, and; 7) Gray Level Dependence Matrix.

Feature selection and RNETS score calculation
Before further analysis, to eliminate the differences in the value scales, we carried out a standardization process for all the extracted image features with z-scores. For feature selection, Spearman's rank correlation analysis was performed and features with correlation coefficient values higher than 0.9 were considered to be redundant and would be filtered out. The LASSO regression was then employed to determine the features with the highest predictive value for the NETS with the penalty parameter tuning conducted by 5-fold cross-validation. The radiomics signature (RNETS) was constructed through a linear combination of the selected features weighted by their corresponding coefficients in the training cohort. The whole workflow of radiomics model is depicted in Figure 2.

Statistical analysis
All statistical analyses and visualization plots were carried out in R studio (version 4.1.0, http://www.r-project.org) and SPSS 22. Correlations between continuous variables were estimated using the Spearman test. Comparisons between two or more group were carried out using either the Wilcox test or the Kruskal-Wallis test. The nomogram and calibration curves were established by using the 'RMS' package. All P values are two-sided and lower than 0.05 are regarded as statistical significance. Detailed workflow of radiomics model (RNETS). received PD-1 inhibitors as second-line treatment who had undergone loco-regional treatments and targeted therapy (Sorafenib or Lenvatinib) prior to immunotherapy. More detailed clinical characteristics of these two Nanfang Hospital cohorts are listed in Table 1.

Development of scoring model based on NETs-related genes
To better apply the NETs-initiation biomarkers to clinical management of HCC, these 69 NETs-related genes were further filtered by LASSO regression analysis to identify six hub genes with lambda.1se values of 0.3077334 ( Figure S1). These six genes (BST1, IL-6, MAPK3, SELP, SELPLG and TLR4) were used to build a NETs-related risk signature called "NETs score", which was calculated for each patients according to the following formula: NETs score = (-0.164×BST1) + (-0.055×IL6) + (1.419×MAPK3) + (-0.518×SELP) + (-0.714×SELPLG) + (-0.574×TLR4). The NETs score (NETS) was a good indicator for the prognosis of HCC patients in the training cohort, and the area under the curve (AUC) values were 0.836, 0.879, and 0.902 for 1, 2, and 3-year overall survival (OS), respectively ( Figure 3A). The same calculation were performed in cohort 2 and cohort 3 as the validation dataset ( Figures 3B, C). The best NETs score cut-off point of 0.8033756 was used to divide patients into high-NETS and low-NETS groups. The low-NETS group had a survival advantage over the high-NETS group in the training cohort (P<0.001) ( Figure 3D). The similar results were also observed in validation datasets (cohort 2, P<0.001; cohort 3, P<0.001) (Figures 3E, F).

Immune microenvironment heterogeneity underlying the NETS
The results from the CIBESORT algorithm showed that immune-activated cells, including plasma cells, activated CD4 memory cells and CD8 T cells were significantly higher infiltration in the low-NETS group ( Figure 4A). The ImmuneScore, StromalScore, and ESTIMATEScore values decreased as the NETs score increased (Figures 4B-D). Compared with the high-NETS group, the low-NETS group exhibited significantly higher expression levels of immune checkpoint genes ( Figure 4E). Moreover, through GESA (Table S2), several immune response related pathways were also activated in the low-NETS group ( Figure S2). In summary, these results suggest that an immunological 'hot' microenvironment exists in the low-NETS group, which tends to respond better to immunotherapy.

Construction of a radiomics biomarker correlated with the NETS
The Spearman correlation analysis was utilized to filter out redundant image features, following which, four radiomics features were eventually selected using the LASSO regression to construct a predictive radiomics biomarker for NETS status (high-or low-NETS) in the training cohort ( Figures 5A, B). The association between the four radiomics features and six NETs-related genes were depicted in Figure 5C. The RNETS was significantly different between high-and low-NETS groups (P<0.001) ( Figure 5D). The AUC value for distinguishing between high-and low-NETS was 0.853 in the training cohort ( Figure 5E).

Prognostic value of the radiomics biomarker
The RNETS cut-off value of 0.8033756 was used to divide patients into high-RNETS and low-RNETS groups. Patients with low-RNETS were associated with significantly longer survival time in the training cohort (P<0.001) ( Figure 6A). The AUC values of ROC curves were 0.813, 0.718, 0.724 for the 1-year, 2-year and 3year OS, respectively ( Figure 6B). The utility of the RNETS was further evaluated in the cohort. Kaplan-Meier curves showed that the low-RNETS group had significantly longer DFS and OS period as compared to those in the high-RNETS group in cohort 4 (P<0.001) (Figures 6C, D). Multivariate Cox regression analysis showed that the RNETS was an independent risk factor for postoperative recurrence. In addition, the BCLC stage and AFP levels were also significantly positively associated with the DFS in cohort 4 ( Table 2). To broaden its clinical application, the variables were then integrated to develop a nomogram model ( Figure 6E). The calibration curve showed there was a good consistency between the predicted value of recurrence and the actual observed value ( Figure 6F). Importantly, in cohort 4, the AUC values of the nomogram for predicting DFS (0.860) were higher than those of RNETS, AFP and BCLC stage ( Figure 6G), indicating that the nomogram has enhanced power in predicting recurrence as compared to RNETS, AFP levels or BCLC stage alone.

Predictive value of the RNETS for anti-PD-1 immunotherapy response
In cohort 5, the RNETS was significantly lower (mean: 0.0254) in the objective response group (CR/PR) than in the SD/PD group (0.2443) ( Figure 7A). When dividing the patients into different RNETS groups with cut-off point of 0.0804, we found that the objective response rate (ORR) was significantly higher in the low-RNETS group than in the high-RNETS group (27.8% vs 10.8%) ( Figure 7B). Moreover, 77.8% of patients in the low-RNETS group had disease control (PR, CR and SD). In the high-RNETS group, only 55.4% of patients had disease control (Table 3). Notably, Kaplan-Meier curve showed that the RNETS was significantly negatively associated with the PFS (P<0.001) ( Figure 7C). These results indicated that the RNETS correlated with the clinical outcomes of anti-PD-1 immunotherapy. Multivariate Cox regression analyses revealed that RNETS, CRP levels, and AFP levels were independent prognostic factors for PFS (Table 4). To establish and validate a risk-scoring model for PFS, we created a nomogram integrating RNETS, AFP levels and CRP levels ( Figure 7D). Calibration curve demonstrated the predicted value of disease progression was in accordance with the actual observed value ( Figure 7E). Our results confirmed that the predictive value of the nomogram outperformed those of the individual RNETS, AFP levels and CRP levels for predicting PFS in patients receiving immunotherapy ( Figure 7F). Subgroup analyses of treatment lines were performed to assess the performance of RNETS in evaluating the response to anti-PD-1 immunotherapy ( Figure S3

Discussion
As only a fraction of HCC patients have shown impressive efficacy to immune checkpoint blockade (30), there is an urgent need to discover robust predictive biomarkers for tracking patient's response to immunotherapy in advanced HCC cases. Unfortunately, there are currently no reliable predictive biomarkers to support clinicians in predicting which patients would respond favorably or unfavorably to immunotherapy. Although PD-L1 and tumor mutation burden (TMB), are two of the most extensively studied predictive biomarkers for cancer immunotherapy (31)(32)(33), their predictive value in gauging patient responses to immunotherapy in HCC has been relatively limited (34,35).
In cancers, the TIME acts as a major role in tumor metastasis, relapse, and resistance to treatment, due to which it has been under intensive investigation. To dissect TIME subtypes may help identify patients who would be likely to respond to immunotherapy. Since neutrophils are an integral part of the TIME, they could function as bridges between the tumor parenchyma and the immune microenvironment (36). Elevated neutrophil-lymphocyte ratio (NLR) in peripheral blood and high infiltration level of TANs in TIME are correlated with poor prognosis in HCC patients (37,38). In addition, NETs, which are unique derivates of neutrophils, have been linked to progression and metastasis in certain solid tumors (10,11,39), and especially in HCC (12, 40). In this work, we have identified a constitutive NETs-related gene signature (NETS) to predict survival of HCC patients. We have noticed that these patients in high-NETS group were associated with poor prognosis. Neutrophils are the most abundant immune cells in peripheral blood and play a fundamental role in inflammatory responses (41); however, their contribution to immune escape in malignancies is still controversial. Recent work has confirmed that TANs play an immunosuppressive role in primary liver cancer and that targeting TANs could be a potentially effective form of immunotherapy for treating liver cancer (42). In addition, the mechanism by which NETs promote immune escape is gradually becoming clearer; for example, it has been shown that NETs act as physical barriers that cover the cancer cells and shield them from immunotherapy (43). Studies have also shown that NETs can suppress T-cell responses to tumors by inducing metabolic and functional exhaustion in these immune cells (44). Our results also demonstrated that the NETS could be used to define TIME subtypes in HCC. Patients in the low-NETS group likely have an immunological 'hot' microenvironment, which may allow them to respond better to immunotherapy. Using ELISA, IHC, and RNA-sequence data, it is possible to assess the abundance of NETs in peripheral blood and tumor tissues. However, these methods are expensive and complicated, which makes their large-scale application impractical. In China, CT examinations are widely used to diagnose liver cancer even in remote areas with underdeveloped medical resources. Therefore, using CT image-based radiological features for prognosticsespecially if they reflect immune-related responses-will be of great clinical and economic value. In studies related to HCC, radiomics was usually used to predict clinical outcomes after surgery or responses to locoregional therapy (22, 23). However, most HCC patients already suffer from an advanced stage of cancer at the time of diagnosis and usually cannot opt for surgical treatment. For such HCC patients, immunotherapy can prolong their survival time and even help them reach a stage where translational surgery becomes a treatment option. In this study, we have developed a radiomics model-based score which we have named the RNETS, as a biomarker for clinicians to use for a quick prognostic indication of the immunological status of a tumor from CT images. In two internal HCC cohorts, the RNETS performed well enough to predict postoperative recurrence of the tumor and response to anti-PD-1 immunotherapy in HCC patients. The present study found that RNETS is associated closely with clinical benefit in patients received immunotherapy as second-line treatment. However, due to the limited number of first-line immunotherapy cohort, studies with a larger sample size should be implemented to further validate that RNETS is capable of predicting immunotherapy benefit for HCC patients, regardless of the therapy line. There still remain some limitations in this study. First, the sample size in the training cohort from the TCIA database is rather small due to a scarcity of cases that have CT image information. Second, the number of NETs-related genes that have been identified is small and may be insufficient for a comprehensive assessment of NETs formation in TIME. Third, since serum NETs expression had been shown to be indirect means for inferring the NETs formation in tumor tissues, we believe that the results in this study will be more profound if we had integrated the information about NETs abundance in peripheral blood into overall analysis. Fourth, this study only covers NETs-related gene signature and tumor heterogeneity is enormous, future projects based on these results should be developed including various key pathway (for example, ferroptosis, hypoxia, epithelial-mesenchymal transition) related prognostic gene signatures and select the optimal radiomics model for implementation in clinical practice. Finally, this study was retrospective, and the results need to be validated using a prospective research framework.
In conclusion, we have developed a radiomics biomarker, the  Predictive value of RNETS for anti-PD-1 immunotherapy response in cohort 5.  RNETS, which links NETs-related gene expression patterns in the TIME to radiomic features obtained from CT images. This scoring model can effectively and noninvasively predict the clinical outcomes and responses to immunotherapy in HCC patients.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement
The studies involving human participants were reviewed and approved by The Ethical Committee of Nanfang Hospital, Southern Medical University. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

Author contributions
HX and YPZ put forward the ideas of this article. HX, QL and YCZ were primarily responsible for conceptualization, methodology, and writing. JH, YS and MJL were responsible for data collection and data curation. JS, ML, MZ and WL were responsible for data revision. YB, YYZ and YPZ revised the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This study was supported by the grants from the National Natural Science Foundation of China (81772923). The funding agency has no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.