This study is based on a multicentric, pathologically re-evaluated series of LNENs, comprising 58 ACs and 111 LCNECs, collected from our previously published studies [14, 28].

For LCNECs, a total of 148 cases were initially identified from the surgical pathology and clinical databases of five Italian institutions (Fondazione IRCCS Istituto Nazionale dei Tumori – INT, Milan; ASST Spedali Civili di Brescia – Brescia; Fondazione IRCCS Ca’ Granda Ospedale Maggiore Policlinico – Milan; Ospedale Policlinico Ospedale San Martino – Genoa; and Humanitas Research Hospital – HRH, Rozzano), between 1975 and 2016 [14]. After central and blinded pathological review, and exclusion of biopsy-only material and cases with only small-cell morphology, 111 LCNECs were retained [14].

For carcinoids, 370 candidate cases were retrieved from the databases of two institutions (Fondazione IRCCS Istituto Nazionale dei Tumori—INT, Milan and ASST Spedali Civili di Brescia– Brescia) between 1988 and 2018. Following the same centralized review and predefined exclusion criteria (lack of surgical resection, biopsy-only specimens, tumors with poorly differentiated components, or uncertain primary origin), 58 ACs were included [28].

In all cases, patients’ charts and tumor morphology were centrally and blindly reviewed by two expert pathologists (C.C. and M.M.) according to the latest WHO classification criteria [14, 28]. Specifically, tumor identification and subtype characterization were performed through the parallel evaluation of at least four consecutive sections from representative formalin-fixed paraffin-embedded (FFPE) tissue blocks, stained with hematoxylin–eosin (H&E), synaptophysin (Syn), chromogranin A (CgA), and Ki-67. Histological assessment included a detailed analysis of tumor architecture, necrosis, mitotic index, and cellular atypia to confirm the initial diagnoses. To determine the suitability of cases for TME immunohistochemical (IHC) analysis, all 111 LCNECs and 58 ACs FFPE blocks were reviewed, and cases with insufficient material for immunoassay were excluded.

The study was performed according to the clinical standards of the 1975 and 1983 Declaration of Helsinki and was approved by the Ethics Committee of Fondazione IRCCS INT (No. INT 171/16).

The IHC study included detection of the following markers: T cells (CD3, CD8), immune suppression (PD-L1, PD-1), HLA (HLA-DR, HLA-I), and fibroblasts (α-SMA) using the antibodies detailed in Supplementary Table 1.

To characterize the immune infiltration and the broader TME landscape, all markers were assessed in the stroma cells within (intra-tumoral) and at a distance (extra-tumoral) from the tumor. PD-L1,HLA-DR, and HLA-I were also evaluated in neoplastic cells. To ensure consistency and minimize interobserver variability, IHC staining results were independently evaluated by two expert pathologists (M.M and C.C.). For all the markers, quantification (both intra-tumoral and extra-tumoral) was based on cell densities in an area of 1mm2 (cells/1 mm2), counted in areas of strongest labelling (“hot spots”). Exceptions included α-SMA, which was evaluated as the percentage of positive cells rather than cell density. Tumor marker evaluation was performed as follows: PD-L1 tumor expression was assessed according to the tumor proportion score (TPS), which categorizes expression into three levels: negative 0%, 1–49%, and ≥ 50% [29]. HLA-DR and HLA-I tumor expression was assessed using the H-score system (range 0–300), which quantifies staining intensity and the proportion of positive tumor cells. The H-score was calculated by assigning reactivity scores to tumor cells as follows: negative (0), slightly positive (1 +), moderately positive (2 +), and strongly positive (3 +). The final score ranged from 0 (all tumor cells negative) to 300 (all tumor cells strongly positive) [30].

Digital slides of CD3 and CD8 IHC-stained sections were acquired using Aperio ScanScope XT® at 40 × magnification for digital pathology investigations. To analyze CD3 and CD8 tumor-infiltrating lymphocytes (TILs), considering their highly heterogeneous distribution, we conducted whole-slide image analysis. CD3 + and CD8 + TILs were quantitatively assessed, as well as their spatial distribution within the tumor microenvironment, through advanced machine-learning-based digital pathology methods.

Specifically, we tested an open-source analytic pipeline for quantifying CD3 and CD8 intratumor and/or stromal infiltrating lymphocytes from digital whole-slide images [31]. In more detail, QuPath was used for image processing [32], and for each CD3 and CD8 virtual slide, tumor areas were manually annotated by two expert pathologists. The annotated areas were then segmented into 1 mm2 tiles, and pathologists manually distinguished the invasive margin (IM) and central tumor (CT) tiles. For each tile, an automated algorithm identified and isolated true positive lymphocyte staining, filtering out artifacts such as speckles and nonspecific background signals (Fig. 1A-B). The validated lymphocyte signals were then saved, and the image was further segmented into super-pixels, which were subsequently classified into tumor or stromal compartments using a random forest classifier trained on data implemented in R and Python. Consequently, comprehensive intratumor (iTIL) and stromal (sTIL) TILs were quantified depending upon their location within the “tumor” or “stroma” super-pixels in the entire slide (ES) and, importantly, in IM and CT (Fig. 1C). For further technical details on the analytical pipeline, please refer to the work of Yoo et al. [31].

A Annotation of tumor area with outer edges of the invasive margin marked in yellow. B “True” positive lymphocyte automatic quantification. C Tumor (red) and stroma (green) identifications with automatic quantification of intra-tumoral and stromal lymphocytes

For each patient, digital pathology analysis of CD3 and CD8 allowed the extraction of several tumor immune microenvironment (TIME) parameters [31]. Specifically, CD3 (iTIL), CD3 (sTIL), CD3 (iTIL/sTIL), CD8 (iTIL), CD8 (sTIL), CD8 (iTIL/sTIL), and tumor-stroma ratio (TSR) were quantified using four summary measures (mean, minimum, median, and maximum) and two heterogeneity metrics (coefficient of variation [CoV] and quartile deviation [QD]). CD3 and CD8 values were calculated at both the IM and CT using the four summary measures. Similarly, five different values (four summary measures plus an overall measure) were derived for the entire tumor area. A total of 207 TIME parameters were initially generated per patient (Supplementary Table 2). To reduce redundancy, we first removed variables that were highly correlated (correlation > 0.99). In addition, in order to identify the most relevant features for distinguishing WHO histology classes, we applied the Boruta feature selection algorithm, a wrapper method based on the random forest classification algorithm. Boruta is specifically designed to retain all features that contribute significantly to the predictive power of a model while eliminating irrelevant ones [33, 34]. The algorithm works by creating shuffled duplicates (shadow features) of the real features and comparing their importance scores (Z-values) obtained from the random forest model. If a real feature consistently achieves a significantly higher Z-value than the highest-ranked shadow feature across multiple independent tests, it is classified as “important” (green zone) and selected. Conversely, features that fail to outperform shadow features are marked as “unimportant” (red zone) and excluded. By applying this method, we identified a robust set of features with meaningful predictive power for WHO histology classes, optimizing model performance while minimizing noise from irrelevant variables. To further investigate the structure of the selected robust digital TIL features, we applied Principal Component Analysis (PCA), a linear dimensionality reduction technique that projects the original high-dimensional data onto a new set of orthogonal components ranked by the amount of variance they explain. Hierarchical clustering was then performed on the PCA-reduced data to identify sample groupings. This unsupervised method constructs a dendrogram based on pairwise Euclidean distances and a complete linkage criterion, progressively merging the most similar clusters.

The relationships between all intra-tumoral and extra-tumoral TME features assessed by pathologists were analyzed using Spearman’s correlation coefficient. Associations between demographic characteristics and clinicopathological features with histology (ACs vs LNECs) or Cluster (1 vs 2 vs 3) were assessed using the Fisher exact test for categorical variables and the Wilcoxon test or Kruskal–Wallis test for continuous variables. Overall survival (OS) was assessed from the date of diagnosis to the date of death for any cause. Disease-free survival (DFS) was assessed from the date of diagnosis to the date of first relapse, death, or last follow-up, whichever occurred first. OS and DFS curves were drawn using the Kaplan–Meier method. The log-rank test was used to assess the survival difference between patient groups. Pearson’s correlation coefficient was used to evaluate the relationship between intra-tumoral and extra-tumoral CD3 and CD8 densities obtained through automated image analysis and the corresponding quantitative counts assessed manually by pathologists. Data analysis was performed using the R environment for statistical computing and graphics (R Foundation, Vienna, Austria- Version 4.3.3). All tests were two-sided, and p-values < 0.05 were considered statistically significant.