Our study was approved by the Ethics Committee of National Cancer Center/Cancer Hospital, Chinese Academy of Medical Sciences (NCC/CHCAMS) (ID Num: NCC2016YJC-08). It was separated into two parts: the clinical cohort and basic research cohort. We enrolled RCC patients with tumor thrombus who were treated surgically from NCC/CHCAMS from January 2000 to December 2019. Patient consent was not required. Patients were followed every 3–6 months for the first 5 years after surgery and then yearly thereafter. All follow-ups were concluded on May 30, 2021. None of these patients received preoperative chemotherapy, radiotherapy, or targeted therapy.

In the clinical cohort, we investigated the prognosis factor for RCC patients with thrombus. The clinical, pathologic, and preoperative laboratory parameters were abstracted from the clinical database. Clinical variables included gender, body mass index (BMI), age at surgery, and paraneoplastic syndrome. Pathologic information included maximum tumor size, tumor laterality, tumor necrosis, Fuhrman grade, sarcomatoid differentiation, perirenal fat invasion, and lymph node metastasis at surgery. Preoperative laboratory parameters included absolute neutrophil, lymphocyte, platelet counts, Hb, LDH, ALT, AST, IgA, IgG, and IgM. NLR was calculated by dividing the absolute neutrophil count by the absolute lymphocyte count. Cutoff values were selected based on our institutional-specific laboratory guidelines, median values, or through related literature review. The clinical-pathological-blood characteristics of the enrolled patients are shown in Supplementary Material, Table S1 . The clinical endpoint in our clinical cohort was defined as CSS.

To further investigate the underlying mechanisms of the prognosis factor, we collected 10 treatment-naive patients of non-metastatic RCC with tumor thrombus from the basic research cohort. Multiple tumor and tumor thrombus regions and peripheral blood mononuclear cell (PBMC) samples were selected for RNA sequencing (RNA-seq). The detailed information of the enrolled patients is shown in Supplementary Material, Table S2 . All of the specimens were verified by histopathology. To early recognize the patients of poor prognosis, we chose recurrence-free survival (RFS) for the outcome endpoint as the time from diagnosis to documented recurrence or metastasis in our basic research cohort. Patients were divided into two groups according to the 3-year RFS. Patients who reached the endpoint were defined as the “High-risk” group. Otherwise, they were separated into the “Low-risk” group.

Public expression data for a further 39 tumor samples from stage III RCC patients with tumor thrombus in TCGA were obtained from UCSC (University of California, Santa Cruz) Xena (https://xenabrowser.net/datapages/). Gene expression data from TCGA datasets were to construct a prognostic gene set for RCC with tumor thrombus.

Statistical analysis was performed using R (version 3.6.0, http://www.r-project.org). Univariable and multivariable Cox proportional hazards regression model were used to estimate hazard ratios (HRs) via the R package “survival.” In addition, the HR, 95% CI, and statistical significance of each key prognostic factor were estimated and demonstrated using a forest plot via the R package “survminer.” The main packages used included survival, survminer, and dplyr.

QC-qualified RNA samples of the preoperation peripheral leukocytes, tumor, and tumor thrombus tissues from the basic research cohort were analyzed using RNA-seq. Strand-specific mRNA libraries were prepared using a NEBNext^® Ultra™ RNA Library Prep Kit (Illumina; cat. no. E7530L). Subsequently, the eligible libraries were sequenced on an HiSeq XTen platform of 150 bp paired-end reads (Illumina; cat. no. FC-501-2521), yielding approximately 20 million (M) reads per sample. Adapter sequences and low-quality reads were cleaned up using Cutadapt and Sickle (http://github.com/najoshi/sickle/) software. Clean reads were quantified against an Ensembl catalog (GRCh37) using Salmon software at the transcript level (24) and aggregated with the R package “tximport” at the gene level; transcripts per million reads (TPM) values were estimated by Salmon gene expression quantification software, and genes with detectable TPM counts over 80% of the samples were maintained. The proportion of immune cells in PBMC was validated by single-cell RNA sequencing (scRNA) information. After the RNA was extracted in PBMC samples, the SMART-based UMI-corrected TCRB libraries were constructed for sequencing, which was detailed in our previous work (25). Data were dealt with the MIGEC, MiXCR, and CellRanger software and the R packages “Seurat”, “immunarch”, and “scRepertoire”.

The gene sets of NETs were derived from two parts. One was a summary gene set of NETs-associated genes, mainly consisting of the ligands and receptors that promote the NETosis, downstream associated signals, and the molecules that illustrated to adhere to the structure of NETs. The other genes came from the neutrophil-associated genes from the weighted correlation network analysis (WGCNA) in the basic research cohort. In total, we converged the 44 overlapping genes in the NETs-associated gene set for signature constructing.

We acquired the marker gene sets for immune cells (26) and classic oncologic pathways (27) from other articles. We carried out ssGSEA to calculate the enrichment score of each term using the R package GSVA” (28). With the above methods, the counts for every immune cell type and the expression for each classic oncologic pathways in each specimen were referred to as the immune-cell score and pathway-activated score. The counts of NETs calculated hereafter stood for as the NETs-score. The signature of NETs in localized RCC with thrombus was constructed from the 44 NETs-associated genes according to the Cox proportional hazard regression model in TCGA cohort, in which each characteristic gene obtained P value <0.05. Every patient could obtain a NETs-score with the signature of NETs.

In this study, WGCNA was restricted to the top 8,000 varying genes according to their standard deviation. Therefore, WGCNA was constructed to find phenotype-related module and hub genes via the “WGCNA” package and R tutorials (29). Finally, highly similar modules and their relationship to phenotype traits (tissue types and prognosis risk) were identified.

By analyzing RNA-seq of tumor and tumor thrombus with an unsupervised clustering algorithm, we separated diverse status on the basis of the infiltration of immune cells and expression of oncological pathway at a much higher resolution. To distinguish the differentially expressed genes (DEGs) in RCC of progression risk group, expression profiles in the “High-risk” tissue were compared with those in “Low-risk” through the “limma” package in R software. Genes with an adjusted P-value < 0.05 and the cutoff value of ∣log₂ fold changes∣> 2 (log₂FC) were defined as DEGs. GSEA was applied to enrich hallmark gene sets downloaded from the Molecular Signatures Database v6.0 (MSigDB) with “GSVA” package in R software. Input genes were listed in descending order according to the log₂FC values. Genes with false discovery rate (FDR) p-value less than 0.05 and nominal p-value less than 0.05 were considered significantly enriched. To investigate Gene Ontology (GO) of a comprehensive set of functionally annotated genes, the R package “clusterProfiler” was used (30), with a cutoff criterion of adjusted p < 0.05.

IHC was performed on archival formalin-fixed paraffin-embedded tumor and tumor thrombus tissue. Tumor thrombus tissues were stained with the following antibodies to identify specific proteins: anti- citrullinated histone H3 (anti-H3Cit) (citrulline R2 + R8 + R17) antibody (Abcam, ab5103; dilution 1:200) and anti-myeloperoxidase (anti-MPO) antibody (Abcam, ab25989; dilution 1:200). NET formation was visually determined as the percentage of neutrophils identified positive for a DAPI and H3Cit signal (31). IHC images were evaluated with image J software.

A total of 214 patients of non-metastatic RCC with tumor thrombus were treated surgically in the period 2000–2019 from NCC/CHCAMS. Clinical-pathological-blood characteristics of these selected patients in our clinical cohort are summarized in Supplementary Material, Table S1 . Of these patients, 203 had available follow-up data and blood count data to calculate preoperative NLR; therefore, this research focused on these patients. The median [interquartile range (IQR)] follow-up was 46.0 months. Median survival time was 127 months (95% CI: 104.2~150.0). The 5-year CSS was 66.3%.

Univariable Cox proportional hazards regression identified the following risk factors: perinephric fat invasion, paraneoplastic syndrome, blood transfusion, BMI, tumor laterality, tumor size, Fuhrman grade 3/4, tumor necrosis, sarcomatoid differentiation, preoperative platelet count, NLR, hemoglobin level, IgG, and IgA ( Table 1 ). A multivariable Cox hazards regression model identified the following independent prognostic factors: NLR ≥ 4 (HR 2.46; 95% CI 1.18–5.1; P = 0.016), Fuhrman grade 3/4 (HR 4.07; 95% CI 1.80–9.2; P < 0.001), and Tumor laterality (left) (HR 2.05; 95% CI 1.08–3.9; P = 0.028) ( Table 1 ).

To further investigate the mechanism of NLR as prognosis factor, we analyzed the local immune infiltration and oncologic features of tumor and tumor thrombus samples from our basic research cohort. Unsupervised clustering on the basic research cohort showed that the tissue samples were predominantly separated into two clusters: low immune cell infiltration (Low-CI) and high immune cell infiltration (High-CI). Patients in the High-CI cluster were featured with rich leukocyte infiltration, such as CD4+ T cells, CD8+ T cells, as well as myeloid cells of macrophages and Neutrophils ( Figure 1A ). The immune characteristic of tumor was consistent with its paired tumor thrombus. In addition, two clusters of unsupervised clustering were consistent with the progression risk group (see the definition in Methods). That is to say, patients in the “High-risk” group characterized as High-CI cluster had the most immune cells highly immersed including T cells, B cells, and Treg cells, compared to the “Low-risk” group ( Figure 1B ). Besides, the total NK cells were with low infiltration in the “High-risk” group ( Figure 1B ).

To further investigate the oncologic pathway-level analyses, we explored simplified signatures of representative genes associated with angiogenesis, cell cycle, antigen presenting machinery, the complement cascade, EMT feature, and metabolism-related pathways, including Fatty Acid Synthesis (FAS)/pentose phosphate and Fatty Acid Oxidation/AMP-activated protein kinase (FAO/AMPK) signaling ( Figure 1A ). Patient tumors in “High-risk” group were primarily characterized as highly proliferative, with enrichment of “cell cycle” pathway related genes ( Figure 1B ). “Low-risk” group also showed increased FAO/AMPK gene signature expression ( Figure 1B ).

To further understand the relationship between neutrophil and prognosis, the infiltration of neutrophil in tumor and tumor thrombus tissues were analyzed. Patients in the “High-risk” group showed an increase in neutrophil infiltration, especially for tumor thrombus ( Figure 1B ). Tumor-derived granulocyte colony stimulating factor (G-CSF) as a NET-inducing factor, and NETs-associated marker genes including Histone family, peptidyl arginine deiminase 4 (PADI4), and matrix metalloproteinase 9 (MMP9), were significantly overexpressed in the tumor thrombus ( Figure 2A ), rather than tumor, of patients with unfavorable prognosis in the “High-risk” group. GSEA revealed that the leukocytes in the peripheral blood of patients with unfavorable prognosis groups were similar to G-CSF-treated PBMCs ( Figure 2B ).

Next, we confirmed that NETosis was morphologically detectable in the tumor thrombus through testing positive for a DAPI, H3Cit, and MPO signal with IHC. NETosis was identified in the “High-risk” group instead of “Low-risk” group ( Figure 2C ). The quantification of IHC staining results indicated that H3Cit and MPO were significantly expressed in the “High-risk” group ( Figure 2D ). These overall modulations suggest that tumor thrombi were divided into NETosis-thrombi in the “High-risk” group and non-NETosis-thrombi in the “Low-risk” group. And NETs induced by tumor cell-derived G-CSF in tumor thrombus have a mechanistic effect of the poor prognosis with RCC.

In WGCNA analysis, the gene co-expression networks were constructed from the basic research cohort. With each module assigned a color, a total of 16 modules were identified and their association with four clinical phenotypes was analyzed, including tissue type (tumor and tumor thrombus) and progression risk (High-risk and Low-risk). The results of the module-trait relationships are presented in Figure 3A , revealing that the black module was found to have the highest association with tumor thrombus tissues with high risk (black module: r = 0.61, p = 0.002). To investigate the underlying function of the genes in black module for the “High-risk_TT” group, we performed GO analysis. In GO analysis, GO terms such as neutrophil activation, neutrophil-mediated immunity, and neutrophil degranulation were enriched in the black module ( Figure 3B ).

As shown in Figure 3C , 486 genes were found in the black module, and 41 genes in neutrophil-characteristic GO-terms in this module. The other was a summary gene set of NETs-associated genes in prior research (see in Methods). As a consequence, the 44 overlapping genes were converged in the NETs-associated gene set for signature constructing ( Figure 3C ). Based on the data of TCGA cohort, a total of 44 genes related to NETs was adapted to Univariable Cox analysis. As a result, a set of five genes—ATAD3B, ADAM8, F3, PSMD13, and SLC11A1—was initially obtained, with a cutoff criterion of adjusted p< 0.05. Subsequently, we got a NETs-score with ssGSEA analysis of these five genes.

Based on the data of TCGA RNA-seq of non-metastatic RCC with thrombus, NETs-score was adapted to the Cox regression analysis, and each variable was given a regression coefficient β (see Table 2 ). The factors of age and NETs-score were candidate risk factors in univariable Cox regression analysis (p < 0.05). Due to Fuhrman grade as an important prognosis factor in RCC, we added it in multivariable Cox analysis with age and NETs-score. As a result, NETs-score (p = 0.021, HR 6.5) and age (p = 0.001, HR 1.1) were independent prognostic factors identified ( Table 2 and Figure 3D ), indicating the NETs-score was a hazard factor for significantly curtailed survival.

According to the TCR repertoire analysis of PBMC via bulk TCR sequencing for the patients in our basic research cohort, the index chao1, a symbol of diversity and richness, was decreased in the “High-risk” group ( Figure 4A ). Each clonotype in the top 25% of two representative patients in the High-risk/Low-risk group was graphically represented ( Figure 4B ), which illustrated that the high-abundance overwhelming clonotypes increased in the “High-risk” group.

Hence, we further investigated the preoperative PBMC transcriptome levels and immune repertoire via scRNA sequencing between two patients in “High-risk” and “Low-risk” groups ( Figure 5A ). Removing the low-quality cells, we acquired single-cell RNA sequencing (scRNA-seq) profiles from 13,532 cells. After normalization of gene expression and principal component analysis (PCA), we used uniform manifold approximation and projection (UMAP) clustering to separate the cells into 16 clusters ( Figure 5B ). These clusters could be assigned to seven known cell lineages through marker genes ( Figure 5C ): Monocytes (2,373 cells, 17.5%, marked with CD14 and FCGR3A), NK cells (1,280 cells, 9.5%, marked with NKG7, GZMH, and GZMA), Dendritic cells (63 cells, 0.47%, marked with FCER1A); macrophages (40 cells, 0.3%, marked with CD68); platelets (110 cells, 0.81%, marked with PPBP); B cells (1,722 cells, 12.72%, marked with MS4A1, CD19, and MZB1); and T cells (7,944 cells, 58.7%, marked with CD4, CD8A). The proportion of each cell lineage varying, obviously, the counts of lymphocytes including naive B cells and CD4 positive central memory T cells (Tcm) decreased while CD4 and CD8 positive effector memory T (Tem) and Treg cells increased in the patient of “High-risk” group ( Figure 5D ). Integrated with the corresponding single-cell immune repertoire data, we confirmed that the increasingly high proportion of clonal expansion (defined as over one T cells shared the same α-β TCR pair) of the CD4 positive Tem cells occurred in the “High-risk” group ( Figure 5E ). Systemic perturbations of lymphocytes existed in the “High-risk” group of non-metastatic RCC with thrombus, featured with TCR diversity decrease and immune function inhibition, which also explains that NLR≥4 was associated with poor prognosis in the clinical cohort.