Downloading transcriptome sequencing data

The scRNA-seq dataset GSE184362, comprising tumor and adjacent normal tissues from 13 THCA patients, was retrieved from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). Additional GEO-derived transcriptome datasets used for meta-analysis are summarized in Table S1.

Bulk RNA-seq data (HTSeq-FPKM) for 512 THCA and 59 normal thyroid samples, together with clinical annotations, were accessed from the TCGA-THCA project (https://portal.gdc.cancer.gov/). To supplement the number of normal tissues, RNA expression profiles from 278 healthy thyroid samples were obtained via the GTEx portal (https://www.gtexportal.org/home/index.html). Since all datasets were publicly available, no additional ethical approval was necessary (Zhao et al. 2020).

scRNA-Seq data processing

Gene expression matrices were generated using Cellranger (10 × Genomics). Cells were filtered in light of the following thresholds: nFeature_RNA between 200 and 5000, nCount_RNA between 1000 and 20,000, and mitochondrial gene content (percent.mt) below 15%. Data quality was tested via unique molecular identifier (UMI) counts and gene correlation analyses. Principal component analysis (PCA) was done on highly variable genes utilizing the Seurat R package (Hao et al. 2021) for dimensionality reduction. The clustering resolution was determined based on the top 10 principal components (PCs) identified in PCA using the UMAP algorithm. The chosen clustering resolution effectively distinguished major cell types without over-partitioning, resulting in 18 initial clusters. The above-mentioned clusters were further annotated into seven different cell types based on known cell marker genes. The JackStraw and ElbowPlot functions guided PC selection, after which clustering was done utilizing a graph-based method and visualized via UMAP. Cell types were annotated in light of established marker genes (Meng et al. 2021). Pseudotime trajectory analysis was done with the help of Monocle2 (v2.22.0; (Qiu et al. 2017)).

Cellular communication analysis in the TME

The CellChat tool (v1.1.3, (Fang et al. 2022)) was used to analyze cellular communication patterns within the TME. First, the normalized gene expression data were input into CellChat for preprocessing, including steps such as identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData, to ensure result robustness. Next, the aggregateNet and computeCommunProb functions were used to calculate the communication probability and information flow intensity of ligand-receptor (L-R) pairs between cell populations, constructing a cellular communication network in the TME. To further explore key ligand-receptor-related genes, the CellPhoneDB tool (v2.1.4, (Efremova et al. 2020)) was also used. The normalized gene expression matrix was input, and preprocessing was performed to identify receptor-ligand pairs between cells. By calculating the significance (P < 0.05) and interaction strength of receptor-ligand molecular pairs, 41 ligand-receptor genes closely related to TME regulation and treatment response were screened, providing a basis for further mechanistic research.

Gene function enrichment analysis

DC-specific marker genes (n = 833) were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the clusterProfiler R package, with significance set at P < 0.05. GO terms were classified into biological processes (BP), cellular components (CC), and molecular functions (MF) (Yu et al. 2012).

Single gene survival analysis

The prognostic relevance of SPP1, IL1B, CD74, IFNGR1, CXCL16, TNFSF13B, AXL, FPR1, IFNGR2, LGALS9, HAVCR2, CD4, TNF, CXCR3, GAS6, CXCR4, and ANXA1 in THCA was assessed utilizing clinical data from the TCGA-THCA cohort. Overall survival (OS) was selected as the endpoint. Cox proportional hazards models were made with the help of the R package survival, and survival curves were generated employing survminer and ggplot2. Group differences in OS between high- and low-expression cohorts were evaluated by the Log-rank test, with P < 0.05 deemed statistical significance (Zhao et al. 2021).

Correlation analysis

Based on the TCGA-THCA dataset, the expression correlation between FPR1 and ANXA1 in tumor tissues was examined employing the R package corrplot. Additionally, correlations between FPR1 and 47 immune checkpoint-related genes (e.g., IDO1, IDO2, LAG3, CTLA4, TNFRSF9, ICOS, CD80, PDCD1LG2, TIGIT, CD70, TNFSF9, ICOSLG, KIR3DL1, CD86, PDCD1, LAIR1, TNFRSF8, TNFSF15, TNFRSF14, CD276, CD40, TNFRSF4, TNFSF14, HHLA2, CD244, CD274, HAVCR2, CD27, BTLA, LGALS9, TMIGD2, CD28, CD48, TNFRSF25, CD40LG, ADORA2 A, VTCN1, CD160, CD44, TNFSF18, TNFRSF18, BTNL2, C10orf54, CD200R1, TNFSF4, CD200, and NRP1) were analyzed (Xu et al. 2022).

Analysis of differential gene expression

The TCGA-THCA dataset and GTEx dataset were utilized to test the differential expression of FPR1 and ANXA1 in THCA tissues. The analysis was completed utilizing the R software package limma, with a threshold set at |logFC|> 1 and P < 0.05 (Ritchie et al. 2015).

Meta-analysis

Meta-analysis was done utilizing the R package meta, with the standard mean difference (SMD) and 95% confidence interval (CI) as effect size measures. Heterogeneity across studies was evaluated employing the Q-test and chi-square test, with I2 and P values guiding model selection: a fixed-effects model was made when P > 0.05 and I2 < 50%, indicating low heterogeneity, otherwise a random-effects model was utilized.

Subgroup analyses were conducted based on detection methods (RNA-seq vs. microarray). Sensitivity analysis was done utilizing a leave-one-out approach, sequentially excluding each study to assess the robustness of the pooled effect and to identify studies significantly impacting the association between FPR1, ANXA1, and THCA outcomes (Balduzzi et al. 2019).

TME analysis

The TME was tested utilizing the R package estimate by calculating the StromalScore, ImmuneScore, and ESTIMATEScore for each sample. Differences in these scores between high- and low-expression groups of FPR1 and ANXA1 were subsequently compared (Huang et al. 2020).

Analysis of immune cell infiltration

The enrichment scores of 24 immune cell-related functional marker genes were calculated using the R software package GSVA for single sample gene set enrichment analysis (ssGSEA) based on a study by Bindea et al. (Hänzelmann et al. 2013). The markers for the 24 immune cells provided in the Immunity article by Bindea et al. were employed to assess the immune infiltration status in TCGA-THCA data (Bindea et al. 2013).

Analysis of immunotherapy in THCA patients

Immunotherapy data for THCA patients were gained from the TCIA database (https://tcia.at/). The Immunophenoscore (IPS) was utilized to predict responses to cytotoxic T-lymphocyte antigen 4 (CTLA-4) and programmed cell death protein 1 (PD-1) blockade. Differences in IPS between high- and low-expression groups of FPR1 were checked utilizing the R package ggpubr (Zanfardino et al. 2019).

Cell culture

Human PTC cell line TPC-1 (CL-0643, Procell), B-CPAP (CL-0575, Procell), human peripheral blood DCs (CP-H179B, Procell), and human peripheral blood CD3+ T cells (PRI-H-00107, Xqxzbio) were cultured in RPMI 1640 (72,400,120, Gibco) containing 10% FBS (10,437,010, Gibco) and 1% penicillin–streptomycin (15,140,163, Gibco). The cultures were maintained at 37 °C in a 5% CO2 atmosphere. When cell confluence reached 60–80%, passaging and transfection were performed (Chen et al. 2022a, b) (Chen et al. 2022a, b).

Immunoprecipitation

DCs were lysed in RIPA buffer (P0013 C, Beyotime, Shanghai, China), and cellular debris was removed by centrifugation. The supernatant was harvested, and the supernatant of each sample was adjusted to the same concentration. Input controls were reserved; the resultant supernatant was immunoprecipitated with anti-ANXA1 (ab214486, Abcam, 1:2000, UK) or anti-FPR1 antibody (ab113531, Abcam, 1:500, UK) along with protein A/G beads (Santa Cruz Biotechnology) for 2 h. Subsequently, beads were triple-washed, denatured at 100 °C for 5 min, electrophoresed, transferred to nitrocellulose membranes (Millipore, Temecula, CA, USA), and subjected to immunoblotting (Chen et al. 2023).

Preparation of lentiviral vectors

The lentiviral vector pSIH1-H1-copGFP (sh-; SI501 A-1, System Biosciences, USA) was adopted for ANXA1 knockdown, while pCDH-CMV-MCS-EF1α-copGFP (oe-; CD511B-1) facilitated FPR1 overexpression. Lentiviral particles were generated by transfecting HEK-293 T cells (iCell-h237, Cell Applications, Shanghai, China) with a packaging kit (A35684 CN, Invitrogen, USA). After 48 h, viral supernatants (1 × 108 TU/mL) were harvested. Target cells at ~ 40% confluence were transduced with the viral suspension for 8 h, followed by culture in DMEM enriched with 10% fetal bovine serum (FBS). Puromycin (5 μg/mL; A1113803, Thermo Fisher Scientific, China) selection was applied for 4 weeks. Table 1 summarizes lentiviral sequences. Silencing efficiency was verified in DC cells (Figure S1), and the most effective sequence (sh-ANXA1-1) was utilized in subsequent experiments (Vecchi et al. 2018).

Cell Co-culture system

To assess the stimulatory capacity of DC, the transfected DCs from each group were co-cultured with CD3 + T cells at a ratio of 1:10 (DCs: T) in a 96-well plate: 3 × 105 DCs from each group were added along with 3 × 106 CD3 + T cells, and the culture volume was adjusted to 200 μL with RPMI 1640 medium. The cells were then incubated at 37 °C with 5% CO2. After 72 h of cultivation, CD3+ T cell proliferation was tested utilizing the CCK-8 assay. For the evaluation of DCs'phagocytic ability, 1 × 105 cells from each group were suspended in 200 μL RPMI 1640 medium and incubated with FITC-Dextran (100 μg/mL; D1844, Invitrogen™) for 4 h. FITC uptake was quantified via flow cytometry.

To investigate the impact of DC cells on TPC-1 cells, 100 μL (containing 1 × 104 cells) of TPC-1 cells were supplemented to each well of a 96-well plate. Subsequently, CD3 + T cells from the co-culture with respective DC groups were added at a ratio of 1:10 (T: TPC-1) to the TPC-1 cells, and the final volume was adjusted to 200 μL. The co-culture was maintained in a 37 °C, 5% CO2 culture system for 48 h, after which TPC-1 cells were harvested for further analysis (Fang et al. 2021).

Transwell assay

TPC-1 cells were evaluated for migratory and invasive capacities employing 24-well Transwell chambers (8 μm pore size; #3422, Corning, NY, USA). For invasion assays, chambers were pre-coated with Matrigel (#354,277, BD Biosciences, CA, USA), diluted 1:1 with serum-free medium and solidified at 37 °C. Migration assays omitted Matrigel pretreatment.

Cells (2 × 105) suspended in 200 μL serum-free medium were loaded into the upper chambers; 800 μL of 20% FBS medium was placed in the lower wells. After 24 h at 37 °C, membranes were processed by PBS washing, fixation with 4% paraformaldehyde, 0.1% crystal violet staining, imaged in five random fields (Nikon Eclipse Ci, Nikon, Tokyo, Japan), and quantified over three replicates (Hou et al. 2021).

Clonogenic assay

A clonogenic assay was utilized to check the TPC-1 cell proliferative capacity. Cells were inoculated at a density of 2000 cells per 6 cm culture dish. After changing the culture medium, the cells were allowed to grow for 14 days. Subsequently, staining was performed using 0.5% (w/v) crystal violet (C8470, Solebao Technologies, Beijing, China), followed by imaging and quantification of colony formation. Clusters of cells consisting of more than 50 cells were identified as a single clone (Wang et al. 2020a, b, c).

CCK-8 assay

Logarithmically growing cells were seeded at 5 × 104 cells/well in 96-well plates for overnight culturing. Cell viability was examined utilizing the CCK-8 kit (E606335, Sangon Biotech, China) by adding 10 μL of reagent at 0, 24, 48, and 72 h. Following 1 h of incubation at 37 °C in 5% CO2, absorbance at 450 nm was documented with an Epoch microplate reader (BioTek, USA) (Guo et al. 2019; Wang et al. 2020a, b, c; Luo et al. 2018). Each condition was tested in triplicate across three independent experimentations.

Ethical statement

All animal experiments were conducted under protocols approved by the institutional Animal Ethics Committee, in compliance with established guidelines. Every effort has been made to minimize the pain and distress experienced by the animals, as well as to reduce the number of animals required for experimentation. Animal housing, care, and experimentations were done in light of internationally accepted standards for animal welfare. Adequate care has been provided to all animals, and measures have been taken to ensure their proper disposition after the conclusion of the experiments.

Construction of subcutaneous tumor model in Mice

Eighteen 6-week-old male NCG nude mice weighing 20 ± 2 g were procured from Jicuiyaokang. They were housed in SPF conditions with a 12-h light–dark cycle, and a suitable temperature range of 22–25 °C. Prior to cell injection, DCs and T cells, along with TPC-1 cells, were transfected in a 200 μL system containing a ratio of 5 × 105 DCs, 5 × 105 T cells, and 1 × 106 TPC-1 cells. Subsequently, the mice were randomized into the sh-NC + oe-NC, sh-ANXA1 + oe-NC, and sh-ANXA1 + oe-FPR1 groups, with each group receiving a subcutaneous injection of 200 μL of the corresponding cell mixture into the left abdominal wall. Tumor volume was tested every two days. On the 12 th day post-injection, the studied animals were anesthetized with 50 mg/kg pentobarbital sodium (B005, Jiancheng, Nanjing, China), then euthanized by cervical dislocation at the point of neck dislocation, followed by tumor removal, photography, weighing, and storage for subsequent detection of relevant factor expression. To minimize subjective bias, three independent researchers assessed each sample separately (Luo et al. 2023).

RT-qPCR

Tissues were processed for RNA extraction utilizing Trizol (16,096,020, Invitrogen, USA). cDNA was synthesized via RT with the Takara kit (RR047 A, Japan). Gene expression analysis was done via RT-qPCR utilizing TaqMan assays (Applied Biosystems, Foster City, CA, USA), with GAPDH as the endogenous control. Each RT-qPCR experiment was conducted in triplicate wells. Primer sequences are detailed in Table 2. The fold change in target gene expression patterns was computed utilizing the 2−ΔΔCt method (Wang et al. 2023).

Western Blot

Following protein extraction utilizing RIPA buffer, its concentration was checked via BCA assay (20201ES76, Yisheng Biotechnology, Shanghai) and normalized. SDS-PAGE was used for protein separation, followed by transfer to PVDF membranes (IPVH85R, Millipore, Germany) via wet transfer. Membranes were blocked with 5% BSA at ambient temperature for 1 h, followed by overnight incubation at 4 °C with anti-ANXA1 (ab214486, Abcam, UK), anti-FPR1 (PA1-41,398, Invitrogen, USA), anti-BATF3 (ab211304, Abcam, UK), and anti-GAPDH (ab9485, Abcam, UK). After TBST washes, secondary antibody incubation was done with HRP-conjugated goat anti-rabbit IgG (ab6721, Abcam, UK) for 1 h. Chemiluminescent signals were developed and analyzed with ImageJ (NIH) to calculate relative protein expression using GAPDH as the control (Wang et al. 2023). For the full length uncropped original western blots, please refer to Supplemental Material file.

Immunofluorescence staining

After formalin fixation and paraffin embedding, tissue sections underwent deparaffinization, antigen retrieval (10 mM sodium citrate, pH 6.0), and blocking with 10% goat serum. Primary antibodies—anti-ANXA1 (ab214486, Abcam, UK), anti-FPR1 (ab113531, Abcam, UK), and anti-BATF3 (ab302568, Abcam, UK)—were applied overnight at 4 °C. After washing, secondary antibodies (fluorescent anti-rabbit IgG, ab150079, Abcam, UK) were incubated at ambient temperature for 1 h, in the dark. Nuclei were counterstained with DAPI (D1306, ThermoFisher, USA) for 10 min. After final washes, slides were mounted with anti-fade medium and imaged utilizing a fluorescence microscope. Quantification of positive cells was done employing ImageJ or built-in software (Zheng et al. 2022).

Flow cytometry

Tumor tissue was digested in HBSS (with Ca/Mg) enriched with DNase I (10 μg/mL) and Liberase (25 μg/mL) at 37 °C for 30 min with intermittent shaking. The suspension was filtered via a 70 μm mesh, cold PBS-rinsed, and centrifuged at 1000 × g for 5 min at 4 °C. Red blood cells were lysed utilizing ACK buffer (A1049201, Invitrogen), quenched with PBS, and centrifuged again. Cell pellets were resuspended in FACS buffer for downstream analysis. For in vitro-cultured cells, after PBS washing and centrifugation under the same conditions, cells were similarly resuspended in FACS buffer.

The aforementioned cell samples were blocked with goat serum IgG (SL038, Solarbio) for 15 min and subsequently stained with a panel of antibodies according to the experimental requirements. For all channels, the gating of positive and negative cells was based on Fluorescence Minus One (FMO) control. Specific antibody-binding cell populations were identified and quantified by flow cytometry. The antibodies used included BATF3 (ab307471, Abcam), CD1a (ab256268, Abcam), CD86 (ab303578, Abcam), CD4 (ab277931, Abcam), and CD8 (ab237365, Abcam) as reported in (Mokuda et al. 2015).

Statistical analysis

Data were initially summarized utilizing descriptive statistics (mean, median, standard deviation, range). For comparisons between two groups, an independent t-test was applied. A one-way ANOVA was adopted for comparisons across three or more groups, while MANOVA was employed for analyses involving multiple factors. Survival analysis was done employing Kaplan–Meier curves, and differences were examined with the Log-rank test. Pearson's correlation coefficient was adopted for continuous variable relationships. Linear and logistic regression were employed for predictive or exploratory analysis. For non-normally distributed data or unequal variances, non-parametric tests (Mann–Whitney U or Kruskal–Wallis) were applied. All analyses were completed utilizing SPSS (IBM) or R. Statistical significance was made at p < 0.05. Graphs were prepared with GraphPad Prism (GraphPad Software, Inc.).