The human tissue samples from small cell neuroendocrine carcinoma of the endometrium involved in this study were from Clinical and Translational Research Center of Shanghai First Maternity and Infant Hospital, Tongji University School of Medicine, Shanghai, 200120, China. The section samples of small cell lung cancer were from patients who underwent pulmonary surgery at Shanghai Pulmonary Hospital, Tongji University School of Medicine, Shanghai, 200433, China. The human NCI cell lines H446 used in this study was from National Collection of Authenticated Cell Cultures. This study was approved by the Clinical Research Ethics Committee of the College of Medicine, Tongji University.

Tissue samples were collected, fixed and then paraffin-embedded. 5 µm sections were obtained and subjected to histopathology and immunohistochemistry (IHC) analysis. We used Hematoxylin and eosin (HE) staining for histopathological examination and immunohistochemistry for expression analysis, respectively. The “UltraVision Quanto Detection System HRP DAB” IHC Kit (TL-125-QDH; Thermo Fisher Scientific, Waltham, MA, USA) was used for IHC assay according to the manufacturer’s protocol. In this study, the primary antibodies used were as follows: anti-cytokeratin pan (GM351529; Gene Tech, Shanghai, China), anti-cluster of differentiation 56 (CD56) (Kit-0028; MXB Biotechnologies, Fuzhou, China), anti-synaptophysin (Kit-0022; MXB Biotechnologies), anti-chromogranin A (MAB-0707; MXB Biotechnologies), anti-Ki67 (IR626; Dako Products, Santa Clara, CA, USA), anti-vimentin (Y23037; Ventana Medical Systems, Oro Valley, AZ, USA), anti-CD3 (ab16669; Abcam, Cambridge, MA, USA), anti-CD20 (M0755; Dako Products), and anti-programmed death-ligand 1 (ab205921; Abcam). Images were taken using a digital pathology scanner (Pannoramic MIDI II; 3DHISTECH Ltd., Budapest, Hungary).

The single-cell suspension was generated carefully and rapidly to ensure high viability as previously described [11]. Briefly, the tissue was initially subjected to mechanical dissociation and then enzymatic degradation with collagenase type I/II (Thermo Fisher Scientific) and DNAse I (Sigma, St. Louis, MO, USA). Following digestion, the cell suspension was subjected to erythrocytes removing by red blood cell lysis buffer (Solarbio Life Science, Beijing, China) and further live-cell enrichment with Dead Cell Removal Kit (Miltenyi Biotec, Bergisch Gladbach, Germany). Hemocytometer and Countess cell counter were used to count cell number and examine cell viability, respectively. Finally, when the cell viability is bigger than 85%, cell suspension at a concentration of 1 $\times$ 10${}^6$ cells/mL was eligible for subsequent 10$\times$ library preparation and sequencing.

The Chromium Single cell 3’ Reagent v2 Kit (10$\times$ Genomics, Pleasanton, CA, USA) and the Chromium Single Cell Controller Instrument were used to generated 10 $\times$ single-cell library according to the manufacturer’s instructions. After library construction quality control, sequence data were generated using Illumina HiSeq X Ten System (Illumina, San Diego, CA, USA).

We used the Cell Ranger software pipeline (version 3.1.0, 10 $\times$ Genomics, https://www.10xgenomics.com/support/software/cell-ranger) to process the reads, including: mapping reads to the human reference genome (GRCh37) and transcriptome and generating a matrix of gene counts versus cells. Then, we used the Seurat (version 3.1.2, https://satijalab.org/seurat/) R package for downstream analysis. Briefly, we removed low-quality cells by quality control and generated normalized count by “NormalizeData” function. Principal component analysis (PCA) was performed to reduce the dimensionality. Then we used the Seurat “Find Clusters” algorithm and “Run TSNE” function to obtain a graph-based clustering result visualized in 2-dimension. “Find Markers” function was used to find differentially expressed genes (marker genes, p 0.05 and $|$log2foldchange$|$ $>$0.58) in each cell cluster. Cell types were identified using canonical marker genes.

Gene Ontology (GO) enrichment analysis for differentially expressed genes (DEGs) was performed using R package based on the hypergeometric distribution. The z scores were computed from normalized –log10 (p value) generated from the Fisher exact test. GO enrichment graph was generated using R package pheatmap (version 1.0.12, https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap).

We performed dimensional reduction and clustering again to generate subclusters for original cluster 2 and cluster 9. Then, 8 subclusters were obtained, subclusters 1–5 were annotated as neuroendocrine tumor cells. The Monocle2 R package (version 2.9.0) [12] was used to perform the trajectory analysis on the tumor cells, so as to infer the relative differentiation time of each cell based on gene expression. Monocle2 R package functions were used to select ordering genes, reduce dimensional, and finally plot genes in pseudotime. The differentiation trajectory of tumor cells or the evolution of tumor cell subclusters during development can be deduced by pseudotime analysis.

Initial copy number variations (CNVs) for each region were estimated using the inferCNV R package [13]. The CNV of total cell types was calculated by the expression level from single-cell sequencing data for each cell with a cutoff of 0.1. Genes were sorted based on their chromosomal location, and a moving average of gene expression was calculated using a window size of 101 genes. The expression was then centered to zero by subtracting the mean. The neuroepithelial (cluster 2) cells were selected as malignant cells, leaving all remaining cells as the normal cells. De-noising was carried out to generate the final CNV profiles. The relative CNV value on each chromosome in each cell was showed in the form of a heat map, and 5 CNV groups was clustered according to CNV level.

TCGA Uterine Corpus Endometrial Carcinoma (UCEC) patient RNA-seq data were downloaded from the TCGA database project by UCSC Xena (https://xenabrowser.net/datapages/). The top and bottom 10% of patients based were selected on ISL1 mRNA abundance followed with differential gene analysis with R package DESeq2 (version 1.42.0, https://bioconductor.org/packages/release/bioc/html/DESeq2.html).

Spearman’s correlation analysis was performed to assess the relationship between subclusters of single cell RNA-seq data and cancer cell line encyclopedia (CCLE) H446 RNA-seq data.

The human NCI cell lines H446 which we used in the study was maintained in our laboratory from Chinese Academy of Sciences. H446 were validated by STR profiling and tested negative for mycoplasma. Then H446 cells were cultured in RPMI Medium 1640 with L-Glutamine (HyClone, Chicago, IL, USA) supplemented with 10% fetal bovine serum and 1% Pen/Strep (100 U/mL penicillin/100 µg/mL streptomycin; Gibco, Carlsbad, CA, USA) under the conditions of 37 °C and 5% CO${}_2$.

H446 cells were seeded in 96-well plates (3 $\times$ 10${}^3$ cells/well) for 24 h. The cell numbers were determined by using Cell Counting Kit-8 (CCK-8) at 450 nm with a microplate spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). As for migration, H446 cells (8 $\times$ 10${}^4$ cells) were seeded into the upper chambers. After incubation for 18 h, each chamber was stained by using AM (calcein-acetoxymethyl ester, calcein-AM; Thermo Fisher Scientific) [14]. At last, the migrated cells were imaged by fluorescence analysis (Nikon, Tokyo, Japan).

After cell counting, the cells were incubated with ISL LIM Homeobox 1 (ISL1) monoclonal antibody (1:200, 15661-1-AP; Proteintech, Rosemont, IL, USA) at 4 °C in the dark for 60 min, followed by incubation for 30 min with a fluorescent secondary antibody (1:200, SA00013-2, CoraLite488-conjugated goat anti-rabbit IgG (H+L)). The cells were washed twice with PBS at 4 °C (1500 rpm, 5 min) for testing. Finally, the positive cells were sterile sorted using the BD FACSAria high-speed flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA).

Data were analyzed with GraphPad Prism (Version 9.5.1, https://www.graphpad.com/) software. For comparisons between two groups, statistical evaluation was done using the two-tailed Student’s t-test. For multiple comparisons, two-way Analysis of Variance (ANOVA) test was used. The association between ISL1 subtype and clinical characteristics were explored by chi-square test. For all statistical tests, p 0.05 was considered statistically significant. Error bars show the standard error of the mean.

A 60-year-old woman presented to our hospital with abnormal vaginal bleeding for more than 1 month. Pathological examination of endometrial curettage showed endometrial malignancy. Therefore, the patient underwent surgery, and pathological examination revealed high-grade SCNEC of the endometrium (Fig. 1A). IHC analysis of the sample showed as following: synaptophysin${}^{++}$, CD56${}^{+}$, chromogranin A${}^{+}$, and Ki-67 80%, immune checkpoint marker Programmed death-ligand 1${}^{+}$ (PD-L1${}^{+}$). CD3 and CD20 staining showed immune infiltration in the tumor environment (Fig. 1B–F). Then we extracted tumor tissue for further single-cell sequencing.

Fig. 1.

Pathology details of the patient with endometrial small cell neuroendocrine carcinoma. (A) Diagnosis and treatment of the patient with endometrial small cell neuroendocrine carcinoma. (B) Hematoxylin and eosin (HE) staining and Pan-cytokeratin (CKpan) expression detected by immunohistochemistry on formalin-fixed paraffin-embedded slides. Scale bars, 80 µm. (C) Expression of phenotype-related genes (CD56, Synaptophysin, Chromogranin A). Scale bars, 60 µm. (D–F) Expression of proliferation gene (Ki67), stromal cell markers (Vimentin, CD3, CD20) and immune checkpoint marker (Programmed death-ligand 1, PD-L1) detected by immunohistochemical analysis. Scale bars, 80 µm.

The primary tumor sample was obtained after resection and subjected to single-cell RNA library preparation using the 10$\times$ Genomics Chromium platform (Fig. 2A). The mean reads per cell was 63,620, and the median number of genes detected per cell was 1248 (Fig. 2B). After quality filtering using the Seurat package, 10,214 cells remained for downstream analysis (Fig. 2C). Graph-based clustering revealed 12 clusters, and known marker genes were used to identify seven major cell types: neuroendocrine tumor cells (neurogenic differentiation 1 [NEUROD1] and neuronatin [NNAT]; cluster 2) [15], fibroblasts (type III collagen [COL3A1]; cluster 9), T cells (CD3 delta subunit of the T-cell receptor complex; clusters 1, 3 and 8), natural killer cells (granulysin [GNLY]; cluster 7), B cells (CD79A; clusters 5, 6, and 10), myeloid cells (CD68; clusters 4 and 11), and platelets [16] (triggering receptor expressed on myeloid cells-like [TREML1]; cluster 12) (Fig. 2D–F, Fig. 3E and Supplementary Fig. 1). IHC staining was performed to confirm the existence of stromal cell types such as fibroblasts, T cells, and B cells with vimentin, CD3, and CD20 antibodies (Fig. 1E). Overall, this analysis showed a complex cellular ecosystem in the SCNEC tissue.

Fig. 2.

Single cell RNA sequencing (ScRNA-seq) profiling of tumor and stromal cells from the small cell endometrial neuroendocrine carcinoma (SCNEC) sample. (A) Experimental workflow of scRNA-seq procedure for the SCNEC tumor. (B) Quality control and metric information for the sequencing sample. (C) The remaining cells after quality control and filtering step. (D) The t-distributed stochastic neighbor embedding (t-SNE) projection where cells that share similar transcriptome profiles are grouped by colors representing unsupervised clustering results. (E) The t-SNE plot demonstrates the major cell types. (F) Expression of representative marker genes of the stromal cell types.

Fig. 3.

Diversity of stromal cell subtypes and expression profile of tumor cells. (A) t-SNE plot of 1451 tumor cells and 8763 stromal cells, color-coded by their associated cluster or the assigned subtype. (B–E) t-SNE plot, color-coded for relative expression (lowest expression to highest expression, gray to red) of marker genes for the T-cell (B), myeloid (C) and B-cell (D) subtypes as indicated. (E) Expression of representative marker genes of the tumor cells. MKI67, Antigen Kiel 67; MZB1, Marginal zone B; MS4A1, Membrane Spanning 4-Domains A1; NNAT, Neuronatin; NEUROD1, Neurogenic differentiation 1; NES, nestin; SNAP25, Synaptosome Associated Protein 25; NCAM1, Neural Cell Adhesion Molecule 1; SYP, Synaptophysin; CHGB, Chromogranin B; CHGA, Chromogranin A; ISL1, ISL LIM Homeobox 1; PAX6, Paired Box 6.

To decode the immune composition complexity, we further identified the cell subsets in lymphocytes and myeloid cells. We annotated two T cell subsets including conventional CD4${}^{+}$ T cells (CD4; cluster 1) and CD8${}^{+}$ T cells (CD8B; cluster 3) (Fig. 3A,B). Proliferative immune cells (marker of proliferation Ki-67 [MKI67]) consisted of both T and B cell lineages (Fig. 3B). We identified two subsets (clusters 4 and 11) of myeloid cells. Dendritic cells (cluster 11) and macrophages (cluster 4) were characterized by the enriched expression of CD1E and CD163, respectively (Fig. 3C). For B cells, follicular B cells (membrane spanning 4-domains A1 [MS4A1]; cluster 10) and plasma B cells (marginal zone B1 [MZB1]; clusters 5 and 6) were further revealed (Fig. 3D). Then we focused on the expression profile of the tumor cells. Differential gene expression analysis of each cluster revealed that cells in cluster 2 exhibited relatively high expression of neuroendocrine relevant genes, including NNAT, NEUROD1, nestin (NES), synaptophysin, and neural cell adhesion molecule 1 (Fig. 3E). These results generate a baseline tumor and stromal cell diversity atlas based on the single-cell NEC transcriptomes.

We identified tumor cell subclusters to explore their functions. The 1835 cells from clusters 2 and 9 were re-clustered into eight distinct subclusters (Fig. 4A). Subclusters 6, 7, and 8 corresponded to endothelial cells (decorin [DCN]), fibroblasts (COL3A1), and epithelial cells (keratin 8 [KRT8] and epithelial cellular adhesion molecule [EPCAM]), respectively (Fig. 4A,B). Subclusters 1–5 referred to cancer cells (Fig. 4B,D). Proliferative markers such as MKI67 were abundant in subcluster 1. Neuroepithelial markers were abundant in subclusters 1, 3, and 5 (Fig. 4B and Supplementary Fig. 2). Subcluster 5 showed the high expression of genes related to nervous system development such as NES and synaptosomal-associated protein, 25 kDa (SNAP25) (Fig. 4B). By contrast, cells in subclusters 2 and 4 showed very limited expression of the neuroepithelial markers (Fig. 4B and Supplementary Fig. 2). The Monocle 2 algorithm was employed to characterize cancer cells. The pseudotime trajectory showed that subclusters 3 and 5 exhibited similar gene expression patterns (Fig. 4C). Subclusters 3 and 5 were both enriched in ISL1 (Fig. 4B, Supplementary Table 1). Studies have shown that abnormal expression of ISL1 is closely related to the occurrence and progression of various cancers such as gastric and prostate cancers [17, 18].

Fig. 4.

ScRNA-seq further analysis identifies distinct populations of tumor cells. (A) tSNE plot of tumor cells and fibroblasts, color-coded by their associated cluster (left) or the assigned cell types (right). (B) Violin plots displaying the expression profile of representative known markers recently reported. (C) Pseudo-time analysis of tumor cells inferred by Monocle2. Each point corresponds to a single cell, and each color represents a tumor subcluster as indicated. (D) Heatmap showing large-scale copy number variations (CNVs) for individual cells (rows). (E) Violin plots showing distributions of CNV scores among different cell clusters and tumor subclusters.

Furthermore, we inferred large-scale chromosomal CNVs in each single cell based on the average expression patterns across intervals of the genome. It is worth noting that the tumor cells have a higher CNV compared with other cell types (Supplementary Fig. 3, Supplementary Table 2). The CNV landscape distinguished malignant cells in subgroups/subclusters 1, 2, 3, 4, and 5 (Fig. 4D). Compared to normal cells, according to the cell chromosome fragment variation, the malignant cells were clustered into five groups, corresponding to the five subclusters (Fig. 4D). We found that cancer cells from subclusters 3 and 5 exhibited higher CNV levels than those from subclusters 2 and 4 (Fig. 4E), suggesting that ISL1-positive cells represent a more aggressive state. Therefore, we carefully investigated the enriched pathway based on ISL1-positive population makers and closely examined the proliferation and migration properties of the ISL1-positive population in the subsequent experiments.

We performed gene ontology enrichment analyses for tumor cell subclusters. Up-regulated genes expressed in ISL1-positive population (subclusters 3 and 5) were mainly enriched for nervous system development-related pathways, such as olfactory nerve development and synaptic membrane (Fig. 5A). Up-regulated genes expressed in subcluster 1 were mainly related to cell cycle. Up-regulated genes expressed in subclusters 2 and 4 were mainly enriched in ribosome functions. Furthermore, we extended the pathway analysis to TCGA UCEC RNA-seq dataset. We selected the top and bottom 10% of patients based on ISL1 mRNA abundance and performed differential gene analysis with DESeq2. This analysis identified 829 upregulated and 1648 downregulated genes between ISL1-low and ISL1-high groups. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis identified ‘Neuroactive ligand-receptor interaction’ as the top enriched term (Fig. 5B–D, Supplementary Table 3). In contrast, Motor proteins related terms are enriched in the downregulated genes (Fig. 5C, Supplementary Table 3). These results suggested that ISL1 might regulate the expression of neuroactive-related genes.

Fig. 5.

ISL1 modulates neuroactive pathway in subclusters of single cell RNA-seq data and The Cancer Genome Atlas (TCGA) patients. (A) The enriched gene ontology terms for marker genes in each tumor subcluster. (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of upregulated genes in ISL1-high compared to ISL1-low patients. (C) KEGG analysis of downregulated genes in ISL1-high compared to ISL1-low patients. (D) Heatmap shows relative abundance of 53 genes in the KEGG term “Neuroactive ligand-receptor interaction” in ISL1-high and ISL1-low patients from the TCGA The Cancer Genome Atlas Uterine Corpus Endometrial Carcinoma (UCEC) cohort.

Large-scale, multicentric, multi-omics analyses of multiple types of cancer provide evidence that SCNC of different tissue origins have similar characteristics [19]. Because there is no mature neuroendocrine tumor cell line in the field of endometrial cancer, in order to further study the effects of ISL1 on the function of NET cells, we used the classic neuroendocrine tumor cell line H446 [20]. H446 is an established human small cell lung cancer cell line. Small cell lung cancer is an extremely aggressive neuroendocrine tumor. Correlation analysis showed that there’s a high correlation between H446 gene expression profile and the neuroendocrine endometrial single cell RNA-seq data (Supplementary Fig. 4). We sorted the H446 cells by flow cytometry according to the expression level of ISL1 on the surface of the tumor cells, and obtained two groups of cells with high and low ISL1 expression for further experiments (Fig. 6A). First, we verified the proliferation function of the two groups of cells. We concluded from the CCK-8 assay that the proliferation ability of cells with high ISL1 expression was higher than that of the low ISL1 expression group (****p 0.0001; Fig. 6B). The results of the EdU cell proliferation assay confirmed these results (***p 0.001; Fig. 6C,D). We assessed the migration ability of the cells, and found that the proliferation ability of cells with high ISL1 expression was significantly enhanced, and the number of migrated cells quantified by the transwell experiment was much higher than that of the low ISL1 expression group (***p 0.001; Fig. 6E,F).

Fig. 6.

Flow cytometric sorting of neuroendocrine tumor cells to obtain ISL1 expression group, and explore the differences in cell function in vitro. (A) Flow cytometric sorting of neuroendocrine tumor cell line H446 according to the difference in the expression of ISL1 on the cell surface. (B) The proliferation ability of H446 cells with high expression of ISL1 is stronger than that of low expression group, (measured by cell counting kit-8 assay, ****p 0.0001, t test). (C,D) The proliferation ability of H446 cells with high expression of ISL1 is stronger than that of low expression group (verified by EdU experiment, ***p 0.001, t test). Scale bars, 200 µm. (E,F) Transwell experiments showed that the cell migration ability of the high expression ISL1 group was stronger than that of the low expression group (***p 0.001, t test). Scale bars, 100 µm. (G) Representative immunohistochemistry (IHC) images of ISL1 in small cell lung cancer (SCLC, n = 30) clinical samples. Scale bar = 5 µm. (H,I) The relation between the expression level of ISL1 and TNM stage, *p 0.05, chi-square test. EdU, 5-Ethynyl-2’-deoxyuridine.

Since small cell neuroendocrine endometrial carcinoma is very rare, we thus further collected 30 SCLC patient tumor tissues and classified into early stage (stage I, II) and advanced stage (stage III, IV) according to National Comprehensive Cancer Network (NCCN) stage. IHC analysis displayed various ISL1 expression levels in these SCLC patient tissues (Fig. 6G). In addition, we correlated the ISL1 expression level with the clinicopathological characteristics of SCLC patients via a chi-square test. Our results showed that the expression level of ISL1 is correlated with SCLC stages (Fig. 6H,I). All the above experiments and analysis demonstrate that tumor cells with high ISL1 expression have stronger malignant behavior than the low expression group, suggesting that ISL1 may have great research value and significance in NETs.