Abstract

   The tumor ecosystem of papillary thyroid carcinoma (PTC) is poorly
   characterized. Using single-cell RNA sequencing, we profile
   transcriptomes of 158,577 cells from 11 patients’ paratumors,
   localized/advanced tumors, initially-treated/recurrent lymph nodes and
   radioactive iodine (RAI)-refractory distant metastases, covering
   comprehensive clinical courses of PTC. Our data identifies a
   “cancer-primed” premalignant thyrocyte population with normal
   morphology but altered transcriptomes. Along the developmental
   trajectory, we also discover three phenotypes of malignant thyrocytes
   (follicular-like, partial-epithelial-mesenchymal-transition-like,
   dedifferentiation-like), whose composition shapes bulk molecular
   subtypes, tumor characteristics and RAI responses. Furthermore, we
   uncover a distinct BRAF-like-B subtype with predominant
   dedifferentiation-like thyrocytes, enriched cancer-associated
   fibroblasts, worse prognosis and promising prospect of immunotherapy.
   Moreover, potential vascular-immune crosstalk in PTC provides
   theoretical basis for combined anti-angiogenic and immunotherapy.
   Together, our findings provide insight into the PTC ecosystem that
   suggests potential prognostic and therapeutic implications.

   Subject terms: RNA sequencing, Cancer genomics, Thyroid cancer, Tumour
   heterogeneity
     __________________________________________________________________

   The characterisation of the papillary thyroid carcinoma (PTC) tumour
   microenvironment remains crucial. Here, the authors perform single-cell
   RNA sequencing in 11 patients and identify potential opportunities for
   the use of immunotherapy and its combination with anti-angiogenic
   therapy in PTC.

Introduction

   The incidence of thyroid cancer has increased by 3% annually in the
   United States over the last four decades, driven largely by the rise in
   papillary thyroid cancer (PTC)^[64]1. Although most PTCs present an
   indolent clinical course, some of them have developed to locoregional
   or even distant metastatic disease at diagnosis. In the recurrent or
   metastatic settings, combination of surgery, radioactive iodine (RAI)
   ablation and thyroid stimulating hormone (TSH) suppression can still
   achieve a favorable prognosis for most cases, while a fraction of
   patients would ultimately progress into RAI-refractory (RAIR) status or
   even succumb to this disease, who may be potential candidates for
   alternative treatments including molecular targeted inhibitors and
   immunotherapies^[65]2,[66]3. The evolving trends of progressive PTCs
   have challenged the clinical practice and promoted researchers to
   further decipher their biological architectures.

   Bulk sequencing of PTC have advanced our understanding of its genetic
   characteristics^[67]2,[68]4,[69]5. For instance, detection of
   BRAF^V600E and TERT promoter mutations can help distinguish malignant
   thyroid nodules and identify patients with dedifferentiation
   potential^[70]5,[71]6. However, the biological underpinnings of PTC
   evolution from early to advanced stage, or from RAI-avid to RAIR state
   remain unclarified. Although bulk sequencing can delineate the genetic
   landscape of the whole tumor entity, it inevitably averages the
   expression profiles of diverse cells and masks the critical differences
   between tumor components. This highlights a critical need to elucidate
   the compositions, properties and underlying mechanisms in the complex
   tumor microenvironment (TME). Single-cell RNA sequencing (scRNA-seq),
   which enables us to quantify features of individual cells, is a
   powerful tool for the investigation of the cellular components and
   their interactions in the TME. Currently, scRNA-seq has been widely
   applied in a broad spectrum of cancers^[72]7. However, improved
   characterization of PTCs at single-cell resolution is still lacking.

   In this work, we perform scRNA-seq to analyze 158,577 cells from 11 PTC
   patients’ paratumors, localized or advanced tumors, initially-treated
   or recurrent lymph nodes (LNs), and RAIR distant metastases, covering
   comprehensive clinical courses of this disease, filling the current
   blank of single-cell profiling of human PTC. Using this unique
   resource, we analyze cell lineages, transcriptional states,
   developmental trajectories and cell-cell crosstalk in PTCs, thereby
   shedding light on the tumor ecosystems underlying PTC initiation and
   progression.

Results

A single-cell expression atlas of papillary thyroid cancer ecosystems

   To comprehensively resolve the tumor ecosystem heterogeneity during PTC
   initiation and progression, we used scRNA-seq (10X Genomics) to profile
   tumor and stromal cells of 23 fresh samples from 11 patients
   (Fig. [73]1a), including six paratumor tissues, seven primary tumors,
   eight involved LNs, and two RAIR subcutaneous loci belonging to distant
   metastases in PTC. In addition to classical PTCs, three patients were
   diagnosed with follicular variant (FV, Case 5 and 7) or tall-cell
   variant (TCV, Case 6) after careful pathological review. Detailed
   clinicopathological information is provided in Supplementary
   Data [74]1. Moreover, the genomic mutations for these patients were
   also assessed by whole-exome sequencing (WES) and Sanger sequencing,
   and the status of key PTC-related mutations (BRAF, RAS, TERT promoter)
   was provided (Supplementary Data [75]2; Supplementary Table [76]1).
   Except for Case 4, 6, and 7 with only recurrent loci, the remaining
   eight patients had paired samples collected for scRNA-seq analysis. For
   example, for case 11, the involved LNs and subcutaneous distant loci
   were identified by computed tomography (CT) scan, clinical examination
   and pathological review to assure the correct tissues for analysis
   (Supplementary Fig. [77]1a, b). By this means, our cohort covered a
   relatively comprehensive collection of tissues mirroring tumor
   progression process, including paratumors, primary lesions, nodal
   metastases and RAIR distant metastases.

Fig. 1. Expression profiling of 158,577 single cells in PTCs.

   [78]Fig. 1
   [79]Open in a new tab

   a Workflow of sample composition, processing and bioinformatic analyses
   for 23 samples in the present study. b t-SNE plot of all high-quality
   cells profiled in the present study colored by major cell lineage. c,
   Heatmap of the canonical and curated marker genes for major cell
   lineages. d t-SNE projection showing the landscape of immune cells,
   colored by cluster (left) and tissue (right). e Tissue preference for
   each immune cell subcluster estimated by Ro/e. Source data are provided
   in the Source Data file.

   A total of 158,577 single cells with a median of 1,215 expressed genes
   passed the stringent quality filtering and were incorporated in further
   analysis (Fig. [80]1b; Supplementary Table [81]2). After integrating
   the transcriptional data from all acquired cells, we primarily applied
   low-resolution t-distributed stochastic neighbor embedding (t-SNE)
   clustering and identified six main cell populations, which were labeled
   as T/natural killer (NK) cells (CD3D, CD3E, CD3G, CD247), B cells
   (CD79A, CD79B, IGHM, IGHD), thyrocytes (TG, EPCAM, KRT18, KRT19),
   myeloid cells (LYZ, S100A8, S100A9, CD14), fibroblasts (COL1A1, COL1A2,
   COL3A1, ACTA2) and endothelial cells (PECAM1, CD34, CDH5, VWF)
   (Fig. [82]1c; Supplementary Fig. [83]1c). Each of these populations was
   captured from different tissue types of different patients
   (Supplementary Fig. [84]1d, [85]e). In addition, all these cell types
   were further validated using another scRNA-seq study of PTC
   (Supplementary Fig. [86]1f, [87]g)^[88]8.

   Subsequently, we aimed to depict a more detailed immune landscape
   through a high-resolution t-SNE analysis. The T, NK, B, and myeloid
   cell lineages were demarcated into 22 finer subclassifications based on
   their patterns of differentially expressed genes (DEGs) (Fig. [89]1d).
   Specifically, T cells were dichotomized according to the expression
   level of CD4 and CD8, which were further divided into 7 clusters for
   CD4^+ cells and 4 clusters for CD8^+ cells, while the NK, B, and
   myeloid cells were divided into 2, 3, and 6 subsets, respectively
   (Fig. [90]1d). Each cluster was assigned with a putative identity
   demonstrating its potential functional capabilities and demonstrated
   different tissue enrichment preferences, as quantified by the ratio of