A B cell-IgA-epithelial axis enhances antitumor immunity and improves outcome in HPV-associated penile squamous cell carcinoma

Domain: neuroscience
Published: 2025-12-11
Editorial outcome: major_revision

Abstract
Penile squamous cell carcinoma (PSCC) is a rare, aggressive malignancy, with HPV driving approximately 50% of cases and correlating with improved outcomes. In this study, we perform single-cell RNA sequencing on tumor microenvironment (TME) cells from 23 treatment-naive PSCC patients, including 11 cases reported previously, and validate our findings in a larger cohort of 85 patients. Our analysis reveals significant immune cell infiltration in PSCC, with HPV+ tumors exhibiting elevated levels of tumor-infiltrating B cells, plasma cells, tertiary lymphoid structures, and IgA secretion. A subset of malignant epithelial cells in HPV+ PSCC expresses the polymeric IgA receptor (PIGR), facilitating IgA transcytosis and inducing anti-tumor responses. Functional studies show that PIGR overexpression promotes tumor apoptosis and enhances immune responses, with these effects diminish in IgA-deficient mice. Elevated IgA and PIGR expression correlate with improved survival. These findings provide a single-cell atlas of the PSCC TME, highlighting the B cell–IgA–PIGR axis as a pivotal driver of antitumor immunity and clinical outcomes in HPV+ PSCC.

Introduction
Penile cancer is a rare urologic carcinoma with orphan disease designation and a prevalence of 0.1–1 per 100,000 men in high-income countries1. However, the incidence is higher up to 10% among male malignancies in the developing countries of African, Asian and South American2. There are various risk factors that cause the different incidence in global regions, including the absence of protective circumcision, chronic inflammatory conditions and HPV infection3,4. The majority of penile cancers originate from the squamous cells and are collectively termed penile squamous cell carcinomas (PSCC), with only about half of patients surviving beyond 5 years1,5. Although rare, PSCC is closely linked to HPV infection, with a pooled HPV prevalence of 50.8%4,6,7. Clinically, PSCC can be divided into HPV-driven and non-HPV-related subtypes8. HPV exerts its oncogenic influence through the viral oncoproteins E6 and E7, which disrupt cell cycle regulation and promote genomic instability9,10. Despite these oncogenic mechanisms, HPV+ PSCC is associated with better outcomes compared to HPV- cases11–13. The reasons for this paradox, however, remain limited.

To date, multi-omics approaches have shed light on the genomic, epigenetic, and transcriptomic landscape of PSCC, revealing somatic mutations14–17, genome-wide methylation patterns, and gene expression profiles18 in PSCC. Recent whole-exome sequencing in PSCC has shown that HPV- PSCC tends to be associated with TP53 loss-of-function mutation16. Another work described the single-cell composition of the PSCC TME and unveiled TP53 mutations as a driver of an aggressive phenotype, irrespective of HPV status19. While these studies have focused on tumor-intrinsic factors, the influence of the immune microenvironment, particularly in HPV+ PSCC, remains unclear.

HPV infection creates a distinct inflammatory immune environment, which may reshape the cellular composition and functional dynamics within the PSCC TME. Given the rarity of PSCC and the limited availability of clinical samples, a detailed characterization of the TME, specifically comparing HPV+ and HPV- PSCC, is essential. Such insights could help clarify the differential outcomes between these two groups20. Furthermore, the validation of these findings across larger PSCC cohorts is needed to ensure robustness and clinical relevance.

In this study, we apply scRNA-seq to explore the heterogeneity and differences of the TME in HPV+ versus HPV- PSCC. We aim to identify key immune cell populations and their roles within the TME. We employ a range of validation techniques, including immunohistochemistry (IHC), multiplex immunofluorescence (mIF), flow cytometry (FCM) and enzyme-linked immunosorbent assay (ELISA) to investigate the immune cell composition and immunoglobulin (Ig) production across a larger PSCC cohort. Our study also explores the heterogeneity of malignant epithelial cells, focusing on a subset expressing the polymeric IgA receptor (PIGR). Using PIGR-overexpressing PSCC primary cells and in vivo xenograft models, we investigate the roles of an immune-PIGR axis in mediating anti-tumor effects.

Results
Single-cell RNA-seq profiling of PSCC
To explore the complexity and heterogeneity of the tumor microenvironment (TME) in PSCC, we performed scRNA-seq on 23 tumor specimens and 6 paired normal adjacent tissues (NATs) from patients as the discovery cohort, including 11 cases reported in our previously published work20 (Fig. 1A, Supplementary Fig. 1A and Supplementary Table 1). One NAT sample (P4N) failed library quality control, leaving 28 high-quality samples (23 tumors and 5 NATs) for downstream analysis. HPV status for the discovery cohort was assessed using p16 immunohistochemistry (IHC)19,21, HPV16/18 DNA in situ hybridization, and HPV PCR (Supplementary Fig. 1B–D). In addition, we collected a larger cohort of PSCC cases to support validation using orthogonal assays (see “Methods”), with detailed clinical and pathological information, as well as sample allocation and time-sequential accrual (Supplementary Fig. 1A and Supplementary Table 2, 3).

Among the 28 samples analyzed by single-cell RNA sequencing, a total of 12.3 billion mapped reads were generated, with 29,845 genes detected (unique molecular identifiers (UMIs) > 0). On average, each sample produced over 441 million sequencing reads with an average sequencing saturation (covering the fraction of library complexity) of 63.2% (Supplementary Table 4). A total of 202,424 single cells passed quality control filters (based on gene expression and mitochondrial genes ratios; see “Methods”) (Supplementary Fig. 2A). On average, each cell had 2209 detected genes and 7138 UMIs, indicating sufficient transcriptome coverage and a representative capture of expressed genes (Supplementary Table 5).

Single-cell landscape of the tumor microenvironment in PSCC
We next used Seurat22 R package for unsupervised clustering and the Uniform Manifold Approximation and Projection (UMAP) for dimensionality reduction and cell visualization (see Methods). We analyzed the total 202,424 cells from the discovery cohort, which captures the full diversity of TME cells, including both immune and non-immune populations, across tumor and NAT samples. All cells expressed high levels of housekeeping genes, such as ACTB, B2M and GAPDH (Supplementary Fig. 2B) and were then grouped into 17 distinct clusters using unsupervised Louvain clustering. Based on known cell markers, we classified 17 clusters into four major cell types: immune cells (PTPRC), epithelial cells (KRT5), endothelial cells (VWF), and cancer-associated fibroblasts (CAFs) (COL6A2) (Fig. 1B, C, Supplementary Fig. 2C). These cell types were consistently identified across PSCC patients, though their relative proportions varied, likely due to patient-specific factors, inherent sampling variability, etc (Fig. 1D). Among them, the identified immune cells emerged as the predominant population within tumor tissues, indicating significant immune infiltration in PSCC (Fig. 1E). Consistent the results from scRNA-seq analysis, mIF confirmed the distribution of immune cells (CD45), epithelial cells (KRT5), endothelial cells (CD31), and CAFs (VIM) across NATs and tumor tissues in PSCC (Fig. 1F).

To validate our scRNA-seq findings, we performed mIF staining on the same PSCC samples, targeting five immune subpopulations, including CD4 and CD8 for TIL-Ts, CD20 for TIL-Bs, CD86 for antigen-presenting cells (APC, macrophages and dendritic cell), and CD56 for NK (Fig. 1G and Supplementary Table 3). In comparison to NATs, the absolute numbers of CD4+, CD8+, CD20+, CD86+, and CD56+ cells per tissue area were significantly increased in tumor tissues (Fig. 1H). These results support our conclusion that PSCC tumors exhibit enhanced immune cell infiltration compared to NATs.

To contextualize our findings, we compared our PSCC data with publicly available datasets from other squamous cell carcinomas (SCCs), including head and neck squamous cell carcinoma (HNSCC)23, nasopharyngeal carcinoma (NPC)24 and cutaneous squamous cell carcinoma (cSCC)25. Similarly, immune (PTPRC), epithelial (KRT5), endothelial (PECAM1) and CAFs (COL6A2) were identified using established markers in each dataset (Supplementary Fig. 2D). While the four cell types were present in all squamous cell malignancies examined, the relative proportions varied across different SCC types (Fig. 1I). Notably, both PSCC and NPC exhibited higher immune cell infiltration, with PSCC showing a dominant immune cell compartment (74.62%) within TME (Fig. 1J).

Distinctive TIL-Bs composition in HPV+ vs. HPV- PSCC
To explore the differences of tumor infiltrating immune cells between HPV- and HPV+ PSCC, we focus on 12 immune clusters from 23 PSCC patients, including 87,790 cells from HPV- tumor and 46,994 cells from HPV+ tumor (Fig. 2A). The remaining immune cell subpopulations were categorized into tumor-infiltrating B cells (TIL-Bs), tumor-infiltrating T cells (TIL-Ts), natural killer (NK) cells, and myeloid-derived cells (Fig. 2A). Among the immune cell populations, we identified two major TIL-Bs populations: SDC1+ plasma cells (C1-Plasma-SDC1) and CD20+ B cells (C2-B-MS4A1). TIL-Ts were subdivided into five clusters (C3-C7), including CD8+ TIL-Ts (C3-CD8T-CD8A), two subtypes of CD4+ TIL-Ts cells (C4-CD4T-FOXP3 and C5-CD4T-CD40LG), a proliferative TIL-Ts cluster (C6-T-MKI67), and early T-cell precursors (C7-PreT-PTCRA). Myeloid-derived cells were grouped into four clusters, including C9-Myeloid-SPP1, neutrophils (C10-Neu-HCAR3), C11-Myeloid-LYZ and mast cells (C12-Mast-MS4A2) (Fig. 2A). For 12 clusters, we found that SDC1+ plasma cells (C1) and CD20+ B cells (C2) were increased in HPV+ PSCC, while cluster 6 was enriched in HPV- PSCC (Fig. 2B, C). Among four immune cell subpopulations, TIL-Bs also emerged as a significant increasement in HPV+ compared to HPV- tumor (Fig. 2D). Interestingly, this enrichment of TIL-Bs has also been observed in other HPV-associated cancers, such as HPV+ HNSCC26–28 and HPV+ oropharyngeal squamous cell carcinoma (OPC)29, suggesting a common feature of HPV-associated squamous cell carcinoma.

To validate our scRNA-seq findings, we performed FCM analysis on CD45+ immune cells isolated from available fresh tumor tissues of 9 HPV+ and 11 HPV- PSCC patients (Fig. 2E, Supplementary Fig. 1A and Supplementary Table 3). The results demonstrated a significant increase in the proportion of TIL-Bs within the total CD45+ immune cells in HPV+ PSCC compared to HPV- PSCC (Fig. 2F). In contrast, the proportions of total CD45+ immune cells, CD3+ TIL-Ts, CD4+ TIL-Ts and CD8+ TIL-Ts showed no significant differences between the two groups (Fig. 2F). We further validated these findings through double-plex IHC staining for B cells (CD20) and T cells (CD3) in tumor sections from our validation cohort of 85 PSCC patients (33 HPV+ and 52 HPV- cases) (Supplementary Fig. 1A and Supplementary Table 3). Importantly, the clinical outcomes observed in our validation cohort were consistent with previous studies, showing that HPV+ PSCC patients generally have better survival than HPV- patients11–13 (Fig. 2G). This consistency in clinical outcomes reinforces the relevance of using this validation cohort for further analysis.

Tertiary lymphoid structures (TLSs), which are immune cell aggregates of B cells and T cells that form outside secondary lymphoid organs (SLOs) in response to chronic inflammation or infection30,31, have been associated with improved survival outcomes in several cancers, including melanoma32, renal cell carcinoma33, and sarcoma34. In our validation cohort, TLSs were significantly more enriched in tumor regions compared to adjacent non-tumor areas (Fig. 2H, I and Supplementary Fig. 3A). Notably, the density of TLSs per tumor area was significantly higher in HPV+ PSCC compared to HPV- cases (Fig. 2J). Further analysis showed a significant enrichment of CD20+ B cells within TLS in HPV+ PSCC tumors compared to HPV- cases, while CD3+ T cell numbers within TLSs remained similar between the two groups (Fig. 2K). mIF staining for CD3, CD20, and CD38 confirmed that TLSs in HPV+ PSCC contained significantly higher numbers of CD20+ B cells and CD38+ activated B cells/plasma cells compared to TLSs in HPV- PSCC (normalized by tumor area to account for sample variability) (Fig. 2L, M). These findings indicate that HPV+ PSCC is associated with a higher prevalence of TIL-Bs and TLSs, contributing to a distinct immune microenvironment compared to HPV- cases.

IgA levels are increased in HPV+ PSCC and correlate with improved survival
To further characterize TIL-Bs in PSCC, we bioinformatically isolated 12,724 TIL-Bs from 23 tumor samples, including 4459 cells from HPV- tumor and 8,625 from HPV+ tumor (Fig. 3A). Based on the expression of key marker genes, we identified four distinct TIL-B clusters in PSCC: BACH2+ or BACH2- Memory B cells cells (C1-C2, CD79A+MS4A1+SDC1-CD24+), plasma cells (C3, CD79A+MS4A1-CD24-CD38+SDC1+) and plasmablasts (C4, CD79A+MS4A1-MZB1+CD38+MKI67+) (Fig. 3A–C). Given the observed higher proportion and number of TIL-Bs in HPV+ PSCC (Figs. 2D, 3D), we next investigated the functional impact of this increased B cell infiltration in HPV+ PSCC. To this end, we performed differential gene expression analysis of TIL-Bs between HPV+ and HPV- PSCC from the scRNA-seq datasets (Fig. 3E). TIL-Bs were divided into HPV+ and HPV- groups to minimize analytical bias. This analysis revealed significantly higher expression of immunoglobulin-related genes such as IGHG2, IGLC2 and IGHA1 in TIL-Bs of HPV+ PSCC compared to HPV- PSCC (Fig. 3E). The top differentially expressed genes were primarily associated with biological processes related to immunoglobulin-mediated and humoral immune responses (Fig. 3F).

Recent studies have demonstrated that plasma cells in tumors produce immunoglobulins that shape anti-tumor immune responses35, with different immunoglobulin isotypes driving distinct immune responses within the tumor microenvironment35,36. To further investigate this in PSCC, we performed bulk B cell receptor sequencing (bulk BCR-seq) on 20 fresh tumor tissues (P11-P30). Library construction failed in 4 cases due to the limited residual tumor tissues, resulting in 16 samples (8 HPV- and 8 HPV+) available for downstream analysis. Following quality control and data processing (see Methods), we found that IgG (70.44%) was the dominant immunoglobulin isotype, followed by IgA (25.66%), while IgE, IgD, and IgM were negligible (Supplementary Fig. 4A). Notably, the diversity of IgG, IgA, IgE and IgM IGH isotype was significantly higher in HPV+ PSCC compared to HPV- PSCC (Supplementary Fig. 4B).

To validate these findings, we performed IHC staining for SDC1 (CD138, a plasma marker), IgA and IgG in our validation cohort of 85 PSCC patients (Supplementary Figs. 1A, 5A, C and Supplementary Table 3). IHC results were scored on a scale from 0 to 4, with higher scores reflecting higher staining intensity and the fraction of positive cells. The detailed scoring system used for CD138, IgA, and IgG in our study was summarized in (Methods). Notably, the proportion of CD138+ signals (score ≥3) was significantly higher in HPV+ PSCC (28/33; 84.84%) compared to HPV- PSCC (14/52; 26.92%) (Fig. 3G). Similarly, IgA+ signals (score ≥3) was also significantly elevated in HPV+ PSCC (26/33; 78.78%) relative to HPV- PSCC (2/52; 3.84%) (Fig. 3H). In contrast, the proportion of IgG+ signals (score ≥3) was comparable between HPV+ (19/33; 57.57%) and HPV- PSCC (29/52; 55.76%) (Fig. 3I). Besides, the mean IHC score of CD138 and IgA was significantly higher in HPV+ PSCC compared to HPV- PSCC (Fig. 3J, K), while IgG levels remained comparable between the two groups (Fig. 3L). In addition, ELISA was performed on freshly harvested tumor tissues and paired plasma samples from 13 PSCC patients (Supplementary Table 3). This analysis revealed a significant increase in tumor-associated IgA levels in HPV+ PSCC, whereas IgG and plasma IgA/IgG levels showed no significant differences between HPV+ and HPV- patients (Fig. 3M, N). When stratifying the 85 PSCC patients into high (score ≥3) and low (score ≤2) expression groups based on IHC staining for IgA, IgG, and CD138, only IgA expression was positively correlated with improved patient survival (Fig. 3O–Q). These findings suggest that HPV+ PSCC exhibits significantly elevated IgA levels compared to HPV- PSCC, and that elevated IgA is associated with better survival outcomes in PSCC patients.

Distinctive functional phenotypes of epithelial cells in HPV+ and HPV- PSCC
Considering that transformed epithelial cells are the origin of PSCC and the primary targets of HPV infection1,6,7, we next sought to explore the intratumor heterogeneity within the epithelial compartments of PSCC. Our primary objective was to differentiate malignant epithelial cells from non-malignant epithelial ones in PSCC TME. Since malignant epithelial cells typically undergo large-scale chromosomal alterations during tumor progression, we used inferCNV23,37,38 to estimate the sample-wise copy number variations (CNVs) of 12,748 epithelial cells isolated from the scRNA-seq datasets of PSCC tumor tissues and NATs from the discovery cohort (Fig. 1B). Epithelial cells from NATs, which showed no significant CNVs, were used as a reference. Epithelial cells exhibiting high CNVs scores were classified as malignant (Fig. 4A), resulting in the identification of 9255 malignant cells for further analysis (see “Methods”).

Unsupervised clustering of the malignant cells revealed four distinct clusters (Fig. 4B). Among these, HPV+ PSCC exhibited an increased proportion of cluster 5 compared to HPV- group (Fig. 4C). Go analysis of cluster 5 marker genes showed the various biological process, including “regutation of viral lifecycle”,”epithelial cell apoptotic process” and “cell chemotaxis” (Fig. 4D, F). Notablely, PIGR was specifically high expressed in the C5 cluster (Fig. 4D, E). PIGR is a receptor found on the surface of epithelial cells that regulates the transcytosis of dimeric IgA and pentameric IgM, produced by local plasma cells39. Clinically, PIGR expression has been shown to vary across cancer types40–43 and is associated with distinct patient outcomes and response to chemotherapy44. These findings highlight the unique biological behavior of PIGR+ malignant cells (C5) and their potential relevance to HPV status in PSCC.

PIGR+ epithelial cells are enriched in HPV+ PSCC
To validate the enrichment of PIGR+ epithelial cells in HPV+ PSCC, we performed IHC staining to assess PIGR protein levels in tumor sections from 33 HPV+ and 52 HPV- PSCC patients (Supplementary Fig. 5D and Supplementary Table 3). The proportion of PIGR+ signals (score ≥3) was significantly higher in HPV+ PSCC (29/33; 87.87%) compared to HPV- PSCC (1/52; 1.92%) (Fig. 4G). In addition, the average PIGR expression level was significantly elevated in HPV+ PSCC compared to HPV- PSCC (Fig. 4H). Using mIF staining, we further examined the colocalization of IgA and PIGR (Supplementary Table 3). The number of IgA-PIGR colocalization events was significantly increased in HPV+ PSCC compared to HPV- controls (Fig. 4I, J), suggesting an increased interaction between IgA and PIGR in HPV+ PSCC.

To investigate the factors driving PIGR expression in HPV+ PSCC, we isolated malignant epithelial cells from a PSCC patient’s tumor lesions, cultured these patient-derived primary cells in conditioned medium and validated their proliferative and tumorigenic characteristics (Fig. 4K). These cells were treated with various stimuli, including immunoglobulins (dominated IgA and IgG in PSCC) and proinflammatory cytokines (e.g., IFN-γ, TNF, IL-1), based on prior studies45. Among the treatments, IFN-γ induced the highest increase in PIGR mRNA expression (Fig. 4L), which was further confirmed by Western blot analysis showing significant upregulation of PIGR protein levels in response to IFN-γ treatment (Fig. 4M). Similar findings were observed in HeLa and HT29 cell lines (Supplementary Fig. 6A–D).

Given the strong induction of PIGR by IFN-γ in our earlier experiments, we next sought to determine whether IFN-γ expression was elevated in HPV+ PSCC compared to HPV- PSCC. To address this, we performed IHC staining to assess IFN-γ expression in tumor sections from 28 HPV+ and 28 HPV- PSCC patients (Supplementary Fig. 5E and Supplementary Table 3). These cases were selected from the validation cohort based on the availability of suitable FFPE sections, as some samples had been partially or fully exhausted in prior analyses. The IHC analysis revealed a significantly higher proportion of IFN-γ+ cells (score ≥3) in HPV+ PSCC (75%; 21/28) compared to HPV- PSCC (7.14%; 2/28) (Fig. 4N). In addition, the mean IHC score of IFN-γ was markedly higher in HPV+ PSCC (Fig. 4O). Further correlation analysis indicated a positive relationship between IFN-γ and PIGR expression, as well as between IgA and PIGR, but not between IgG and PIGR, based on IHC staining in PSCC tumor sections (Fig. 5P–R). Collectively, these findings suggest that the enrichment of PIGR+ epithelial cells in HPV+ PSCC is likely driven by elevated levels of IFN-γ within the tumor microenvironment, potentially enhancing IgA-PIGR interactions in HPV+ tumors.

High PIGR expression exerts anti-tumor effects in PSCC primary cells
To explore the functional role of elevated PIGR expression in HPV+ PSCC epithelial cells, we considered evidence from other cancer types where reduced PIGR expression is linked to tumor progression, such as gastric adenocarcinoma46, non-small cell lung carcinoma47, and nasopharyngeal carcinoma43. Given that HPV+ PSCC is associated with better clinical outcomes compared to HPV- cases, we hypothesized that higher PIGR expression in epithelial cells may have an anti-tumor effect in PSCC.

To test this hypothesis, we constructed lentiviral PIGR overexpression vectors and stably infected PSCC primary cells (Supplementary Fig. 7A). Following infection, ZsGreen-positive cells were sorted by flow cytometry to select cells with stable and high PIGR expression (Supplementary Fig. 7B). Successful overexpression of PIGR was confirmed through flow cytometry and Western blot analysis (Fig. 5A, B). We then assessed the impact of PIGR overexpression on cell proliferation and migration. Compared to the control cells, PIGR-overexpressing PSCC primary cells showed reduced proliferation and colony formation (Fig. 5C, D) and decreased migration capacity in transwell assays (Fig. 5E). To further investigate the role of PIGR in vivo, we subcutaneously implanted control- and PIGR-overexpressing PSCC primary cells into nude mice and monitored the tumor growth every other day. As expected, PIGR-overexpressing PSCC cells led to significantly slower tumor growth and lower tumor weight compared to controls (Fig. 5F, G). In addition, clinical correlation analysis in our cohort revealed that higher PIGR expression was associated with improved survival outcomes in PSCC patients (Fig. 5H). Taken together, these findings suggest that elevated PIGR expression in PSCC epithelial cells exerts an anti-tumor effect, both in vitro and in vivo, as demonstrated in immune-incompetent mouse models.

High IgA-PIGR expression promotes anti-tumor activity in PSCC primary cells
Given the elevated levels of TIL-Bs, IgA, and PIGR+ malignant epithelial cells in the HPV+ PSCC tumor microenvironment, we hypothesized that the PIGR-mediated transcytosis of IgA might establish an immune-epithelial axis with potential anti-tumor effects. To test this hypothesis, we conducted RNA sequencing on IgA-stimulated PSCC primary cells, comparing PIGR-overexpressing cells with controls. Principal component analysis (PCA) confirmed the high quality and reproducibility of the sequencing data (Fig. 5I). Gene ontology (GO) analysis of differentially expressed genes revealed that IgA-PIGR groups exhibited significant upregulation of pathways, such as ‘positive regulation of programmed cell death’ (Fig. 5J, K). Conversely, genes related to the ‘mitotic cell cycle’ and ‘DNA replication’ were notably downregulated in the IgA-PIGR group (Fig. 5J, L).

To further validate these findings, we performed gene set enrichment analysis (GSEA). The IgA-PIGR group exhibited significant enrichment in immune-related pathways, including ‘interferon alpha response,’ ‘interferon gamma response,’ and ‘apoptosis’ (Fig. 5M, O–Q). In contrast, the control group displayed enrichment in cell cycle-related pathways such as ‘G2M checkpoint’ and ‘E2F targets’ (Fig. 5N, O–Q). In addition, qPCR analysis confirmed that immune-related genes such as IFIT3 and IFIT2 (IFN-stimulated genes), IL23A (interleukin), and CXCL11 (chemokine) were significantly upregulated in the IgA-treated PIGR group (Fig. 5R). Our results align with findings from studies on endometrial48 and ovarian cancer41 cells, where IgA transcytosis via PIGR also triggers immune responses and promotes apoptosis, suggesting a conserved mechanism across cancer types. Notably, in PSCC cells, we observed specific downregulation of pathways such as ‘Cell cycle,’ ‘G2M checkpoint,’ and ‘E2F targets,’ which were not significantly impacted in endometrial or ovarian cancer. This suggests that the IgA-PIGR axis may have additional inhibitory effects on cell cycle regulation in PSCC, contributing to its anti-tumor activity in this context.

To assess the functional role of the IgA-PIGR axis in regulating cell behavior, we performed cell proliferation, migration, and apoptosis assays. PSCC primary cells, either PIGR-overexpressing or control, were treated with IgA, and the effects on cell proliferation were assessed using CCK-8 assays and colony formation assays. IgA-treated PIGR-overexpressing cells exhibited significantly reduced proliferation and colony formation compared to the IgA-treated control group (Fig. 5S, T). Migration assays using a transwell system demonstrated that PIGR-overexpressing cells had significantly reduced migration capacity compared to controls (Fig. 5U).

We extended these findings by investigating the effects of IgA-PIGR interaction in the A549 lung adenocarcinoma cell line, which naturally expresses high levels of PIGR. The results mirrored those in PSCC cells, showing reduced growth and migration upon IgA treatment in PIGR-expressing cells (Supplementary Fig. 8A–D). Moreover, apoptosis assays using PI/Annexin V staining showed that IgA treatment significantly increased apoptosis in PIGR-overexpressing PSCC primary cells compared to controls (Fig. 5V, W). These findings demonstrate that the IgA-PIGR axis has potent anti-tumor effects in PSCC epithelial cells, impacting key processes such as cell proliferation, migration, and apoptosis. This suggests a critical role for IgA-PIGR interaction in mediating anti-tumor responses in the HPV+ PSCC tumor microenvironment.

High IgA-PIGR expression promotes anti-tumor activity in vivo
To assess the anti-tumor role of the IgA-PIGR axis in vivo, we utilized the MC38 colon carcinoma model, as human PSCC primary cells cannot be used in immune-competent C57BL/6 mice. While the MC38 model has limitations due to its low B-cell and plasma cell infiltration49–51, it remains useful for studying the downstream effects of IgA-PIGR interactions on tumor cells, particularly in an immune-competent setting. In light of this, we hypothesized that the IgA-PIGR axis may also play an anti-tumor role in MC38 tumors.

To test this, we generated PIGR-overexpressing MC38 cells using lentivirus-based PIGR vectors. Cells with stable PIGR expression were selected via ZsGreen-positive sorting through flow cytometry (Supplementary Fig. 7C). Successful PIGR overexpression in MC38 cells was confirmed by both flow cytometry and Western Blot (Supplementary Fig. 7D, E). Next, we subcutaneously inoculated control and PIGR-overexpressing MC38 cells into wild-type C57BL/6 mice and administered anti-CD20 antibodies (aCD20 mAb) at the time of tumor implantation and at 10-day intervals to deplete B cells (Fig. 6A). We monitored tumor growth every other day and observed that depletion of CD20+ B cells from the blood and spleen effectively reduced IgA and IgG levels in tumor tissues (Fig. 6E, F). B cell deletion via anti-CD20 treatment at earlier time-point significantly inhibited MC38 tumor growth, consistent with previous studies showing B cell-deficient mice suppresses MC38 tumors52. In wild-type mice, PIGR-overexpressing MC38 cells exhibited significantly slower tumor growth compared to controls, but this growth inhibition was abolished in mice treated with aCD20 mAb (Fig. 6B–D), suggesting that CD20+ B cells contribute to PIGR-mediated anti-tumor activity in the MC38 model.

To specifically assess whether the IgA-PIGR axis mediates anti-tumor effect, we inoculated control and PIGR-overexpressing MC38 cells into wild-type and IgA-deficient (IgA-/-) mice (Fig. 6G). The absence of IgA alone did not significantly affect MC38 tumor growth (Fig. 6H), indicating that IgA may not be the primary driver of tumor suppression in this model. This aligns with the limited B cell infiltration and restricted isotype-switched immune responses typically seen in MC38 tumors, where IgA is more likely sourced from peripheral sites rather than produced locally by tumor-infiltrating B cells. However, our focus was on the downstream effects of IgA-PIGR interaction. We observed enhanced anti-tumor effects with PIGR overexpression, which were notably reduced in IgA-/- mice, indicating the necessity of IgA for PIGR-mediated tumor suppression (Fig. 6H–J). Flow cytometry analysis further showed increased infiltration of cytotoxic CD8+ T cells in WT mice bearing PIGR-overexpressing tumors, an effect not seen in IgA-/- mice (Fig. 6K–M and Supplementary Fig. 9A). These findings underscore the critical role of an intact IgA-PIGR axis in driving anti-tumor responses in the MC38 model, highlighting IgA’s contribution to PIGR-mediated tumor suppression.

High IgA-PIGR expression correlates with improved survival in PSCC and mucosal-associated carcinoma
To explore the clinical relevance of our findings, we investigated whether the elevated levels of IgA and PIGR in the TME correlate with improved patient survival in PSCC and other cancers. We performed IHC analysis of IgA and PIGR expression in tumor tissues from 85 PSCC patients. The results showed that 27.05% of patients (23/85) exhibited high expression of both IgA and PIGR, while 58.82% (50/85) had low expression of both. In addition, 8.23% (7/85) displayed high IgA but low PIGR expression, and 5.88% (5/85) showed high PIGR but low IgA expression (Fig. 6N). Notably, patients with high expression of both IgA and PIGR demonstrated significantly improved survival outcomes (Fig. 6O).

To extend our findings beyond PSCC, we analyzed IgA and PIGR transcript expression in various cancers using data from The Cancer Genome Atlas (TCGA), including bladder, breast, colon, cervix uteri, esophagus, head and neck, and lung cancers. While no significant correlation with overall survival was observed in cancers such as bladder, breast, esophageal, and lung, mucosal-associated carcinomas—such as cervical squamous cell cervical squamous cell carcinoma (CESC) (Fig. 6P), head and neck squamous cell carcinoma (HNSCC) (Fig. 6Q), and colon adenocarcinoma (COAD) (Fig. 6R)—showed a positive association between elevated IgA and PIGR expression and improved survival. These findings suggest that the IgA-PIGR axis plays a role in favorable clinical outcomes in PSCC, particularly HPV+ cases, and may also contribute to improved survival in other mucosal-associated carcinomas.

Methods
Patient cohorts and sample allocation
A total of 107 penile squamous cell carcinoma (PSCC) cases were included in this study. Specimens were obtained from patients who underwent surgical excision at The First Affiliated Hospital of the University of Science and Technology of China. Written informed consent was obtained from all participants in accordance with the Declaration of Helsinki. Of the 107 cases, 60 archival FFPE samples (P31–P90) were collected prior to project initiation (2011–2021) and were primarily used for validation assays, including immunohistochemistry (IHC). In parallel, 47 fresh tumor cases were prospectively and consecutively collected between 2021 and 2025, serving as the main discovery cohort. Among these, 30 patients (P1–P30) contributed both FFPE and fresh tumor tissues, while an additional 17 patients (P91–P107) were accrued during the revision stage and specifically contributed fresh tumor tissues for scRNA-seq. In total, 28 samples (23 tumors and 5 NATs) passed scRNA-seq quality control, yielding 202,424 single cells for downstream analyses. Among the 23 tumor samples included in the scRNA-seq dataset, 11 samples (P1–P6 and P91–P95) have already been reported in our previously published work20.

Tissue allocation to assays was dictated solely by real-time availability and technical feasibility, without any pre-selection based on clinical or biological features. Fresh samples with sufficient viable tissue were prioritized for scRNA-seq, a tissue-consuming assay that precluded parallel flow cytometry. Flow cytometry was therefore performed on other consecutively accrued fresh tumors that met viability thresholds. For histologic validation, we examined an inclusive FFPE cohort of 90 consecutively accrued patients (P1–P90). To ensure accurate HPV classification, five archival cases with discordant results (p16-positive/HPV PCR-negative) were excluded, leaving 85 evaluable cases for IHC analyses.

HPV status definition
HPV status in PSCC samples was determined using p16 immunohistochemistry (IHC) staining according to the recommendations of the International Society of Urological Pathology (ISUP), which recognizes p16 IHC as a standard method for detecting HPV infection and classifying penile tumors21. p16 positivity was defined as strong, diffuse nuclear and cytoplasmic staining in at least 70% of the tumor cells, based on validated criteria for multiple HPV-related cancers, including PSCC55 and HNSCC26,27. For PCR-based confirmation, HPV-specific universal primers SPF1/GP6++ was used to detect a broad range of HPV types, including HPV11, 16, 18, 31, 33, 34, 35, 39, 51, 52, 53, 56, 58, 61, and 66 following established protocols56 (Supplementary Table 6). DNA was extracted from tumor samples using either fresh tissue (n=47) with the Animal Tissues/Cells Genomic DNA Extraction Kit (D1700, Solarbio), or FFPE samples (n=60) with the TIANamp FFPE DNA Kit (DP331, Tiangen). A PCR product with a ~180bp band indicates HPV positivity, while no visible band indicates HPV negativity. Given the limited sensitivity of PCR in degraded FFPE tissues, HPV status was assessed using both p16 IHC and HPV PCR. Cases were classified as HPV-positive if they were p16-positive and showed an ~180bp band on HPV PCR. Cases were considered HPV-negative if they were p16-negative and showed no PCR band.

HPV DNA in Situ Hybridization Detection
HPV DNA was identified using the HPV detection kit (Triplex International Biosciences, HPV high-risk-20A). Tissue slices were first subjected to microwave treatment for 30 min, and then the slices were subjected to proteinase K digestion for 5 min, followed by probe hybridization (add 20 μl of hybridization solution to each slice and add a coverslip). The slices were placed on a 95 °C heater for denaturation treatment for 5 min, and then were placed in a humidified box and incubated at 37 °C for 16 h. After that, diaminobenzidine (DAB) was used to develop color, and finally hematoxylin counterstaining was used to display the cell nuclei, alcohol dehydration, and neutral gum covered the slices.

Immunohistochemistry (IHC)
IHC was performed retrospectively on patient tumor tissue sections against p16, CD3, CD20, CD138 (SDC1), IgA, IgG, PIGR and IFN-γ. Owing to limited biospecimen availability, IHC staining was performed only once for each tumor sample. Tissues were fixed in 4% paraformaldehyde, followed by dehydration, paraffin embedding and sectioning. Slices were subjected to immunohistochemical staining. Specifically, specific marker-stained images were analyzed digitally with the image analysis software ImageJ to determine the amount or intensity of IHC staining by measuring the area of the IHC-stained cells or structure in the center of a tumor and in the invasive margin of the tumor. The tumor region was delineated by two trained pathologists. Specific markers expression was categorized by tumor cell IHC scoring: high expression (score=4 or 3), low expression (score=2 or 1) or no expression (score=0). The following primary antibodies were used for immunohistochemical staining: mouse-anti-human p16 (MXB-bio, MAB-0673), mouse-anti-human CD3 (ZsBio, ZM-0417), mouse-anti-human-CD20 (ZsBio, ZM-0039), rabbit-anti-human CD138 (ZsBio, ZA-0584), rabbit-anti-human IgA (ZsBio, ZA-0446), rabbit-anti-human IgG (ZsBio, ZA-0448), rabbit-anti-human PIGR (Abram, ab96196), rabbit-anti-human IFN-γ (Abram, ab25101).

The multiplex immunofluorescence staining protocol was performed as follows: First, 3-μm-thick slices from formalin-fixed paraffin-embedded (FFPE) tissues were mounted on highly adhesive glass slices. The slices were then baked at 65 °C overnight, dewaxed using xylene, and rehydrated in a series of ethanol solutions. Following this, heat-induced antigen retrieval was carried out with citrate buffer (Akoya Biosciences, NEL811001KT). After a 10 min blocking step, slices were incubated with primary antibodies at room temperature for 60 min. Unbound antibodies were removed by washing with TBST (0.1% Tween 20), and a fluorescein-conjugated secondary antibody was applied. Finally, nuclei were counterstained with DAPI for 5 min, and the slices were mounted with an antifade medium and stored at 4 °C in the dark. Primary antibodies were mouse-anti-human KRT5 (ZsBio, TA800727), mouse-anti-human CD45 (ZsBio, ZM-0183), mouse-anti-human CD31 (ZsBio, ZM-0044), rabbit-anti-human VIM (Abram, ab92547), mouse-anti-human CD3 (Servicebio, GB12014), mouse anti-human CD4 (ZsBio, ZM-0418), mouse-anti-human CD8 (Servicebio, GB12068), mouse-anti-human CD20 (Servicebio, GB14030), mouse-anti-human CD38 (Servicebio, GB114831), mouse-anti-human CD56 (ZsBio, ZM-0057), rabbit-anti-human CD86 (CST, 91882S), rabbit-anti-human IgA (Abcam, ab124716), rabbit-anti-human PIGR (Abcam, ab96196).

Tertiary lymphoid structures (TLSs) definition and assessment
We used CD3 and CD20 immunohistochemistry (IHC) to identify tertiary lymphoid structures (TLSs) in our 85 FFPE PSCC tumor sections. CD3 and CD20 are well-established markers for T cells and B cells, which are the primary components of TLS. The use of CD3/CD20 IHC staining is a widely accepted method for the initial identification of TLS based on the presence and spatial organization of T and B cells32–34. Only 37 slides contained both tumor and adjacent non-tumor regions within the same tissue section from 85 PSCC cases. To ensure consistency and accuracy in our analysis, our criteria for TLS inclusion: (A) Structured organization: TLS were defined by a clear, organized structure, with B cells and T cells forming distinct zones. (B) Size threshold: only TLS consisting of more than 50 cells were included in the analysis to distinguish them from smaller, less organized immune cell clusters. (C) Normalization by tumor area: the number of TLS was normalized by the area of the tumor sections to provide a consistent measure across samples. This approach aligns with established methods for TLS identification and quantification in tissue sections, which often use similar size thresholds and normalization techniques to ensure robust and reliable data57,58.

IHC Scoring of CD138, IgA, IgG, PIGR and IFNγ
Specific CD138, IgA, IgG, PIGR and IFNγ stained images were analyzed digitally with the image analysis software ImageJ. The stratification of CD138 expression in our study was based on a semi-quantitative scoring system that combines staining intensity with the percentage of positive cells. The scoring system used was as follows: (0) Negative: No immune cells positive for CD138. (1+) Low: 1–5% of immune cells positive for CD138. (2+) Moderate: 6–20% of immune cells positive for CD138 with moderate or strong intensity, or 6–100% of immune cells with weak intensity. (3+) High: 21–75% of immune cells positive for CD138 with moderate or strong intensity. (4+) Very High: 76–100% of immune cells positive for CD138 with moderate or strong intensity. This scoring system is consistent with those used in previous studies that assessed CD138+ plasma cell density in tissue sections59,60. For IgA and IgG, the scoring was based solely on the fraction of positive cells because they are secreted proteins, and the fraction of positive cells is a more direct indicator of the presence of these immunoglobulins in the tissue. The scoring system used was as follows: (0) Negative of immune cells; (1+) 1–5% of immune cells; (2+) 6–20% of immune cells, (3+) 21–75% of immune cells and (4+) 76–100% of immune cells. Similar fraction-based scoring for immunoglobulins has also been applied in studies analyzing immunoglobulin expression in tumor tissues59.

For PIGR and IFNγ IHC scoring, we measured the integrated option density (IOD) and area value of the positive signal in each field of view by ImageJ software, and then calculated the mean density, that is, mean density = IOD/area, which reflects the unit area concentration of the target positive signal. Finally, the average value of the mean density of 5 random areas of each sample was taken as the positive signal intensity value of this sample. The cell staining intensity was scored into 5 levels, (0) Negative: no positive staining, and the positive signal intensity values of the remaining samples were sorted and processed into four categories according to the order ratio of the positive signal intensity values: <25% was scored as 1 score, 26–50% was scored as 2 score, 50–75% was scored as 3 score, and >75% was scored as 4 score.

Flow cytometry
PSCC tumor tissues underwent three to five washes with PBS (1% penicillin-streptomycin, Gibco) were minced into small pieces about 1 to 5 mm2, then were digested in 10 mL of 5 mg/mL collagenase type IV (Solarbio C8160-100mg) in Advanced DMEM/F12 (Gibco) containing 2.5% FBS and 1% P/S for 0.5 h at 37 °C gently shaking. Cell cluster-containing supernatant was collected and centrifuged at 800 × g for 5 min, and the cell pallet was washed three times using Advanced DMEM/F12.

For single-cell suspensions of tumor samples, flow cytometry was performed by staining with Zombie NIR viability dye (BioLegend, 1:500), blocking with human TruStain FcX™ (BioLegend, 1:100), and staining for 30 min at 4 °C with the following anti-human antibodies: CD45 (BioLegend, HI30, 1:200), CD3 (BioLegend, OKT3, 1:200), CD19 (BD BioLegend, HIB19, 1:200), CD4 (BD Biosciences, RPA-T4, 1:200), CD8 (BD Biosciences, RPA-T8, 1:200). For intracellular cytokine staining of tumor-infiltrating T cells in MC38 tumor samples, cells obtained from the digested tumor mass were incubated in a tissue culture incubator for four hours at 37 °C with Cell Stimulation Cocktail (plus protein transport inhibitors, eBioscience™, 1:500). Surface staining was stained with Zombie NIR viability dye (BioLegend, 1:500), blocked with anti-CD16/32 (BioLegend, 1:100) and analyzed by flow cytometry (Beckman colter, CytoFLEX S) after incubation for 30 min at 4 °C with the following anti-mouse antibodies: CD45.2 (BioLegend, 104, 1:200), CD11b (BioLegend, M1/70, 1:200), TCRβ (BioLegend, H57-597, 1:200), CD8α (BioLegend, 53-6.7, 1:200), GZMB (BioLegend, QA16A02, 1:200), IFN-γ (BioLegend, XMG1.2, 1:200) and TNF-α (BioLegend, MP6-XT22, 1:200).

For experiments requiring cell sorting, cells were first incubated with viability factor (Zombie NIR™ Fixable Viability kit) in 1× PBS for 15 min, followed by washing with MACS buffer (1× PBS with 0.5% FBS and 2 mM EDTA) and centrifugation at 800 × g for 5 min and then stained against anti-human CD45 antibody in MACS buffer for 30 min. Cells were then centrifuged again, resuspended in MACS buffer, and sorted using a BD FACSAria™ III cell sorter. Sort cells were collected directly in 1× PBS with 0.04% BSA. For Single-cell RNA-seq analysis, live CD45+ cells from PSCC tumor tissues were sorted by using Fixable Viability Dye Zombie Aqua™ (BioLegend) and anti-human CD45 antibody conjugated to Brilliant Violet 605™ (Biolegend, clone HI30).

Cell lines and PSCC primary cells
MC38, HT29, Hela and A549 cell lines were purchased from the American Type Culture Collection. PSCC primary cells was made in our laboratory. All cell lines were proven to be negative for Mycoplasma contamination. All cells were cultured in a 37 °C incubator with 5% CO2. All cancer cells, except for PSCC primary cells, were cultured in high sugar DMEM supplemented with 10% FBS and 1% penicillin–streptomycin. For the establishment and passage of PSCC primary cells, PSCC tumor samples underwent three to five washes with PBS (1% penicillin-streptomycin, Gibco), were minced into small pieces about 1 to 5 mm2 and then digested in 10 mL of 5 mg/mL collagenase type IV (Solarbio C8160-100mg) in Advanced DMEM/F12 (Gibco) containing 2.5% FBS and 1% P/S for 0.5 h at 37 °C gently shaking. Then, cell cluster-containing supernatant was collected and centrifuged at 800 × g for 5 min. Cell pallet was washed three times using Advanced DMEM/F12 and resuspended in PSCC primary tumor cell culture medium (PRS-PECM-2D, Hefei PreceDo Medical Laboratory Co., Ltd.). Fresh medium was changed every 2 to 3 days. After PSCC primary tumor cell clones were formed, enzymatic digestion was used to remove mixed primary fibroblasts at each cell passage. Repeat the fibroblast removal operation 3-5 times. Finally, cellular immunochemistry was used to clarify the type of squamous cell carcinoma (Guangzhou KingMed Diagnostics Group Co., Ltd.). For cell treatment with various cytokines and IgA/IgG, 1 ng/ml IFN-β, 5 ng/ml IFN-γ, 20 ng/ml TNF-α, 10 ng/ml IL-1β, 0.5 μg/ml IgA, and 0.5 μg/ml IgG were used, and the cells were collected 24 hours after stimulation for further RNA isolation. PSCC primary cells, HT29 and Hela cells were also stimulated with 4 ng/ml and 10 ng/ml IFN-γ to quantify PIGR protein level by western blot.

For the generation of stable PIGR expression cell lines in PSCC primary cells and MC38 cells, lentivirus-based plasmids containing human PIGR and mouse PIGR cDNA were constructed. PSCC primary cells and MC38 cells were seeded in 24-well plates and infected with hPIGR/mPIGR lentivirus for 72 h, with 8 μg/ml polybrene (Sigma, cat. no. TR-1003). After 72 hours, the infection efficacy was evaluated by ZsGReen fluorescence intensity. DAPI (BioLegend, 422801) was stained for dead cells, and stable PIGR expression cells were sorted by FACS using DAPI-negative and ZsGReen-positive fluorescence for downstream experiments.

Mice and xenograft tumor models
C57BL/6 wild-type mice (male, 6 weeks old) and BALB/c-Nude mice (male, 6 weeks old) were purchased from GemPharmatech Co., Ltd and Beijing Vital River Laboratory Animal Technology. IgA−/− mice were purchased from Cyagen Biosciences, and the genetic background was C57BL/6J (male & female, 8 weeks old). All mice were housed in rooms on a standard 12 h light–dark cycle, with a temperature range of 68–72 °F and a humidity range of 30–70%. In all tumor model experiments, 6–12-week-old male mice were used.

For xenograft tumor models, 5×105 MC38 cells (ctrl- or PIGR-overexpressed) or 1×107 primary human PSCC cells (wild-type) cells were subcutaneously injected. Tumor volume was calculated as: 0.5×(L×W2), in which L is length and W is width. In accordance with the ethical guidelines, subcutaneous tumors were allowed to reach a size not exceeding 10% of the animal’s body weight, which constituted the permissible limits established by the ethics board. The maximum permitted tumor burden was not exceeded at any point during the experiment. Tumor tissues were dissected mechanically into single-cell suspensions for flow cytometry or retained for RNA and protein isolation. B cells were depleted from C57BL/6 WT mice by intraocular injection of 200 μg of anti-CD20 antibodies (anti-CD20, BioLegend, SA271G2, Cat# 152116, or 1× PBS) at the time of tumor cell implantation and again at 10 days intervals to deplete B cells. Mice would be killed if the tumor size/burden exceeded 2.0 cm in diameter. For tumor TME analysis, implanted tumor tissues were removed and digested in 5 mL of 1 mg/mL collagenase type IV (Biosharp, BS165-100mg) in RPMI 1640 containing 40 μg/ml DNase I (Biosharp, BS137-10mg) for 0.5 h at 37 °C gently shaking. Cell cluster-containing supernatant was collected and centrifuged at 500 × g for 5 min, and cell pallet was washed three times using MACS buffer for downstream FACS sorting, FCM analysis and other cell assays.

Western blotting
All cell samples were lysed using the RIPA Lysis Buffer (Beyotime, P0013B) with 1% PMSF. The lysates were separated by 10% SDS-Page and transferred onto PVDF membranes. The membranes were then blocked with 5% milk for 1 h. Following blocking, the membranes were incubated overnight at 4 °C with the primary antibodies: anti-human PIGR (Abcam, ab96196), anti-mouse PIGR (Abclonal, A6130), anti-HA monoclonal antibody (Proteintech, 66006-2-Ig) and anti-beta actin monoclonal antibody (Abcam, ab8226). Then, the membranes were washed with 1× TBST and incubated with secondary antibodies at room temperature for 1 h. The protein bands were visualized using the ECL chemiluminescence substrate kit.

Enzyme-Linked Immunosorbent assays (ELISAs)
IgA/IgG levels in patient plasma/PSCC tumor tissues were determined using IgA/IgG Human SimpleStep ELISA Kit (Abcam, ab196263 & ab195215). IgA/IgG levels in mouse plasma/subcutaneous tumor tissues were determined using IgA/IgG Mouse SimpleStep ELISA Kit (Abcam, ab157717 & ab157719). Follow the product instructions strictly: a. Add 50 μL of all samples or standard to appropriate wells. b. Add 50 μL of the Antibody Cocktail to each well. c. Seal the plate and incubate for 40 min at room temperature on a plate shaker set to 400 rpm. d. Wash each well with Wash Buffer. e. Add 100 μL of TMB Development Solution to each well and incubate for 5 min in the dark on a plate shaker set to 400 rpm. f. Add 100 μL of Stop Solution to each well. Shake the plate on a plate shaker for 1 min to mix. Record the OD at 450 nm.

CCK-8 proliferation assays
CCK-8 Cell Proliferation and Cytotoxicity Assay Kit (Solarbio, CA1210) was used to detect cell proliferation. The cells were digested into single-cell suspension, and then the cell suspension was seeded in a 96-well plate (2000 cells/well). Place the culture plate in an incubator to pre-incubate for 8 h (37 °C, 5% CO2). Add 10 μl of CCK reagent to each well. Incubate the culture plate in the incubator for 1 h. Use a microplate reader to measure the absorbance at 450 nm. Testing was performed every 24 h until the end of 96 h.

Colony-formation assays
Cells were cultured in DMEM (10% FBS) for 24 h before seeding and then diluted into a concentration of 1×105/ml. A total of 200 cells (2000 of PSCC primary tumor cell) were seeded in 2 ml of DMEM (10% FBS). Treatments were refreshed on day 0 after seeding. On day 12–14, one image was taken per well (three wells per treatment group). Colony number was calculated using ImageJ.

Transwell migration assays
Cells were starved in FBS-free DMEM for 24 h. Add 600 µl of DMEM containing 20% FBS into the lower chamber of the 24-well plate and use tweezers to place the Transwell chamber into the 24-well plate. The starved cells were digested, resuspended in serum-free DMEM, and 10,000 cells (100 ml) were seeded in the upper chamber. The 24-well plate was cultured at 37 °C and 5% CO2 for 24 h. Take out the Transwell chamber and gently wipe the cells inside the chamber with a cotton swab. Add 4% paraformaldehyde to fix the cells for 30 min. Add 0.1% crystal violet staining solution and stain for 10 min. Randomly take pictures of 3 fields of view in each chamber and use ImageJ to count and average the number of cells for quantitative research.

Generation of single-cell RNA-seq libraries and sequencing
For PSCC tumor (P1-P6 and P91-P107) and normal adjacent tissues (P1-P6), whole single cells were centrifuged for 5 min at 800 × g and were resuspended in PBS with 0.04% BSA. Cells were then counted using the Countess® II Automated Cell Counter (Thermo Fisher) with a final viability of >85%. Single-cell RNA-seq libraries were prepared using the Chromium Single Cell 3ʹ Reagent Kits v3 (10x Genomics). Cell suspensions were loaded onto a Chromium Chip along with the reverse transcription master mix and single-cell gel beads aiming for 2000–10000 single cells per reaction. Following the generation of single-cell gel bead-in-emulsions (GEMs), reverse transcription, purified and amplification were performed. The above single-cell RNA-seq libraries were generated with a unique sample index for each sample and sequenced on a NovaSeq platform (Illumia) to generate 150 bp paired-end reads.

Alignment, quality control of scRNA-seq data
Following sequencing, Raw data (raw reads) of fastq files were assembled from the raw BCL files using Illumina’s bcl2fastq converter and processed through primary quality control. All the downstream analyses were based on the clean data (clean reads) with high quality. Reads from clean data were mapped to the human GRCh38 reference genome by Cell Ranger (10X Genomics, version 6.0.1) pipeline using default parameters to generate gene/barcode matrices.

After alignment using Cell Ranger, a unified gene/barcode matrix for each sample was acquired and read into Seurat (version 4.3.0) for downstream analysis. The 28 single-cell RNA-seq libraries were sequenced to average more than 441 million sequencing reads for each sample, with the average of sequencing saturation (covering the fraction of library complexity) at 63.2%. After creation of the gene/barcode matrices, a cell-level filtering step was performed to remove cells that expressed fewer than 200 genes, more than 5000 genes and those with more than 10% mitochondrial gene expression. Next, genes that were lowly expressed (genes expressed in fewer than 10 cells) and mitochondrial and ribosomal genes across each sample were removed. Finally, 202,424 single cells from a total of 28 libraries were retained for downstream analysis, with 7138 mean UMIs and of 2209 mean genes per cell.

Processing and clustering of scRNA-seq data
Gene expression matrices were normalized for library size per cell and log-transformed for better interpretability using the Seurat function ‘NormalizeData’. For data integration and batch correction, ‘FindIntegrationAnchors’ and ‘IntegrateData’ were applied to the data in Seurat. For previously published HNSCC (GSE103322)23, NPC (GSE150430)24 and cSCC (GSE144236)25 datasets, files downloaded from GEO Accession Viewer were processed, imported into a Seurat-compatible object and integrated by tumor types for downstream analysis. After generating the library-size normalized and integrated gene/barcode matrix, we next sought to reduce the dimensionality of the data for downstream visualization and clustering. Highly variable genes (top 2000) were identified and selected by ‘FindVariableGenes’ based on high dispersion and carried forward into scaling and principal component analysis (PCA) using the Seurat function ‘ScaleData’ and ‘RunPCA’. Louvain clustering was performed on all cells with the ‘FindClusters’ function using the first top PCs and various resolution parameters. And we then used the Seurat function ‘RunUMAP’ with the equal PCs number for visualization.

Identification and quantification of cell types in scRNA-seq data
Identification of various cell clusters was performed using the Seurat function ‘FindAllMarkers’ or FindMarkers’ in conjugation with the differentially expressed canonical markers. The differentially expressed genes with false-discovery rate (FDR) ≤ 0.05 and log2-fold change of average expression ≥ 0.5 were considered significantly different in different cell clusters. For integrated data in multi-tumor types, we annotated four major populations in TME according to the classical gene signature. For malignant cells, we used inferCNV (version 1.10.1) to estimate the sample-wise copy number variations (CNVs) following the standard workflow for epithelial cells with default parameters and selected out 9255 malignant cells with high CNVs for further analyses. We then evaluated the enrichment of cells from a given group type in each cluster by a two-step strategy. Firstly, we summated the cells from total clusters in each sequencing sample, and calculated the proportion of each cluster by dividing sum value of per sample. Secondly, analysis of proportion for sample ratio values by student t test was used to determine if the cell enrichment across clusters was statistically significant differences in cell frequencies.

Bulk B cell receptor sequencing (bulk BCR-seq) and analysis
We collected fresh tumor tissues from 20 PSCC patients (PT11-PT30) and performed bulk BCR-sequencing on the Seqhealth Technology platform to evaluate the types and diversity of immunoglobulin (Ig) isotypes in PSCC. Total RNAs were extracted from cryopreserved snap frozen tumor tissues using TRIzol (Invitrogen™, Cat. No. 15596018). About 2 μg total RNA of each sample were used for BCR sequencing library preparation using KC-DigitalTM Stranded BCR-seq Library Prep Kit following the manufacturer’s instruction. The library products corresponding to 250–500 bp were enriched, quantified and finally sequenced on the DNBSEQ-T7 sequencer. Library construction failed in 4 cases due to the limited residual tumor tissue, leaving 16 samples (8 HPV+ and 8 HPV-) for downstream analysis.

Following sequencing, clean reads were first assigned to clusters based on identical UMI sequences. Within each cluster, reads were pairwise aligned, and those exhibiting >95% sequence identity were further extracted into a new subcluster. Once all subclusters were established, we performed multiple sequence alignment to obtain a consensus sequence for each. These non-redundant consensus sequences served as the input for bulk BCR-seq analysis. Using MiXCR software (version 3.0.3)61, the consensus sequences were aligned to the international ImMunoGeneTics (IMGT) database62 to obtain V, D and J gene fragments, rearrangements and CDR3 sequences. We defined a BCR clone according to the same V(D)J genes, CDR3 nucleotide and amino acid sequences. The diversity of IGH isotype was evaluated by Shannon’s entropy, which have been widely used for assessing the diversity of BCR as described63,64.

Survival analysis using The Cancer Genome Atlas
Molecular and clinical data from TCGA for cervical squamous cell carcinoma (CESC), head and neck squamous cell carcinoma (HNSCC) and colon adenocarcinoma (COAD) were downloaded from Genomic Data Commons Xena Hub (https://gdc.xenahubs.net). A total of 521 patients with matched clinical information and tumor RNA-seq data was used in this study, including CESC (n=144), HNSCC (n=293), and COAD (n=84). We evaluated the expression of IgA based on the average of IGHA1 and IGHA2. For CESC and HNSCC data, the cutoff of high or low expression for PIGR or IgA was the median. For COAD, the lower and upper quartiles of PIGR or IgA were selected for clinical grouping. Kaplan–Meier survival analysis and the log-rank tests were performed to compare overall survival between groups.

PSCC primary cells transcriptome sequencing and analysis
PSCC primary cells in culture in PRS-PECM-2D medium were treated with or without 0.5 μg/ml of natural human IgA or 1× PBS for 24 h Total RNA was isolated with TRNzol Universal total RNA extraction kits (TIANGEN, DP424). About 2 μg total RNA were used for stranded RNA sequencing library preparation using KC-DigitalTM Stranded mRNA Library Prep Kit for Illumina® (Catalog NO. DR08502). The library products corresponding to 200–500 bps were enriched, quantified and finally sequenced on DNBSEQ-T7 sequencer (MGITech Co., Ltd. China) with PE150 model.

Following sequencing, raw sequencing data was first filtered by Trimmomatic (version 0.36). Low-quality reads were discarded, and the reads contaminated with adapter sequences were trimmed. Clean reads were further treated with in-house scripts to eliminate duplication bias introduced in library preparation and sequencing. The de-duplicated consensus sequences were mapped to the hg38 reference genome using star (version 2.5.3a) with default settings. Uniquely mapped reads were counted with in-house scripts at the gene-level using GENCODE V38 transcript annotation for human. Principal component analysis was performed with normalized log-transformed data using the R language function ‘prcomp’. Differential expression genes analysis was performed using DESeq2 (version.1.34.0)65. The significantly expressed genes were selected with foldchange ≥1.5 and adjusted p-value ≤ 0.05.

Signature enrichment analysis
Gene set enrichment analysis (GSEA) was performed by using R package clusterProfiler66 function ‘GSEA’ (version 4.0.5) with gene sets (Hallmarks and KEGG) from the Molecular Signatures Database (msigdbr, version 7.5.1). Gene Ontology (GO) enrichment analysis in this study was performed by using clusterProfiler or the Metascape online website67.

Statistics & reproducibility
Data are presented as means ± SEM. Statistical analyses were performed using GraphPad Prism (v8.4.0) and R (v.4.1.2). A two-tailed paired and unpaired Student t test between two groups were used to determine significance. Chi-Squared test was used to determine significance for the proportion of high (score ≥3) or low (score ≤2) IHC score between HPV+ and HPV- groups. For time-course analysis of tumor growth and CCK8 cell proliferation, a two-way ANOVA test was used. Survival was analyzed with the two-sided log-rank (Mantel–Cox) test. Statistical parameters were represented in the Figure Legend of each Figure. Pearson correlation coefficients was computed to evaluate the correlation between two variables. P<0.05 was considered significant. *P<0.05, **P<0.01, ***P<0.001, and ****P<0.0001. For scRNA-seq data, DEG analysis was performed by using the Wilcoxon rank-sum test. Other statistical analysis used in this paper are described in the corresponding methods. The western blots, HE and fluorescence images shown were selected as representative of results from at least three independent replicate experiments.

Ethics approval
All human samples, plasma, tumor and adjacent tissues, were collected in accordance with the Declaration of Helsinki. The study design was approved by the Ethics Committee of First Affiliated Hospital of the University of Science and Technology of China. The approval ID is 2016KT67. All animal experimental protocols were approved by the Institutional Animal Care and Use Committee of the University of Science and Technology of China. The approval ID is USTCACUC192401055.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "I would like to salute the authors efforts in conducting the study titled “A B cells-IgA-epithelial axis enhances antitumor response and improves outcome in HPV-associated penile squamous cell carcinoma”. I have major clinical and methodological concerns for this study as it would benefit from answering the below major comments:\nAbstract:\n- please clarify : “we collected a cohort of 90 PSCC patients and performed single-cell RNA analyses on intratumoral cell populations from 10 38 treatment-naive primary PSCC cases.” Why only 10 samples from the 90 PSCC pts collected were used any selection bias to clarify?\n- Please avoid the use of words such as “massive” “remarkably enriched” without values to support.\nIntroduction: would benefit from improved referencing and structing such as:\n-avoiding to use concluding statements in the introduction before the methodology and results presentation “We revealed that the intra-tumoral B cells 92 (TIL-Bs), plasma cells, tertiary lymphoid structures (TLSs) and IgA production were 93 4 significantly enriched in HPV+ PSCC. Interestingly, a subset of HPV+ PSCC malignant 94 epithelial cells expressed polymeric IgA receptor (PIGR) to mediate IgA transcytosis, and 95 elicited transcriptional changes that contributes to hindering malignant progression. 96 PIGR-overexpression promotes tumor cell apoptosis and enhances anti-tumor response 97 in PSCC primary cells and in xenograft tumor models, while the phenotype was abolished 98 in IgA deficient mice. In PSCC and some other tumor types that derived from mucosa, 99 the higher IgA and PIGR expression in tumor tissues correlated with improved survival 100 outcomes. Collectively, our work firstly established the comprehensive picture of the 101 PSCC ecosystem, and demonstrated that the B cell-IgA-PIGR axis contributes to the 102 clinical advantage in HPV+ PSCC.”\nMethodology:\n- In the abstract it is stated that they performed scRNA seq on 10 samples in the methodology discovery cohort used is 6 samples only, please correct and clarify\n- The validation cohort of 90 pts they did not perform scRNA seq for validation for results furthermore the second cohort is from the same institution? Need to clarify.\n- Ned to clarify the methodology for calling p16 IHC positive vs negative this is not clearly detailed in the methodology section and is the main message of this manuscript. Multiple quantitative and qualitative methodology for p16 IHC calling can be used and this is a major important limitation for the validation cohort as no HPV PCR nor ISH was performed.\n- The TLS methodology for staining and calling needs the supporting best practice methodology reference.\n-\nResults:\n- Figure 1 needs clarification on the methodology to highlight the number of cases and the location of the tumor collection penile primary alone or with LN and visceral metastatic sites.\n- The number of pts is not clear they state 6 treatment naïve SCC in the results section.\n- Using HPV 16 and 18 DNa testing is not sufficient to detect all possible HPV subtypes, this has potential major implications on the results as the paper is focused on the differences between HPV + and -.\n- Table 1Supplement has major issues:\nMissing HPV status per PCR only presented data on p16IHC, will need to include multiple methodology and to assess the discrepancy between p16 IHC and HPV PCR or other HPV validated testing methods.\n-in table S1 the author have major errors that raise concern on the validity of the data:\nThey have 42 patients with positive tumors on Lymph node dissection so we would expect a minimum of 42 patient with N1 or N2 disease they only include 15 pts in the N stage.\nAlso, they have 17 pts as stage 3 which requires N1 or N2 which would be an error with the 42 positive LND and also with the 8 pts with N1 and the 9 pts with N2.\n- Fig S1 a P1 T3N0M0 is placed as stage III and Pt19 t3n0M0 is stage II\n- Would recommend not splitting into discovery cohort 1 and 2 best to have one discovery cohort and perform similar comprehensive testing. It is not feasible to make definitive conclusions since the discovery cohort 1 only included 2 HPV negative pts.\n- Fig 1 G the used immune markers are not appropriate to evaluate the immune microenvironment, additional markers like Cd8, Cd4, Cd86, Cd263, cd168 and other immune markers of interest would have been more representative of the author’s message.\n- “we bioinformatically isolated 12,314 immune cells from scRNA-seq datasets of discovery cohort 1, including 6,029 cells from HPV- tumor and 6,285 cells from HPV+ 212 tumor (Fig. 2A)” Need to clarify the methodology for this sentence to select the immune cells as this cohort only has 2 hpv negative vs 4 HPV + tumor samples, this raises concerns of selection bias and quality of data analysis performed thereafter in that section of the manuscript.\n- From the discovery cohort 2 of 90 pts the authors need to highligtht the reason to select only 9 HPV + and HPV negative SC to validate in fig 2 D-E\n- Clarify selection of the N=37 for fig 2G\n- Fig 2 H can be moved to supp\n- Fig 2 I and J have errors and conflicting results with text please clarify.\n- Please explain the selection of only 4 HPV = and 4 HPV – pts for bulk BCR seq, none of the cohort has 4 HPV – and 4 HPV + pts?\n- Please add reference supporting the 5 levels for Cd138 stratification\n- Please add references and explain the different scoring system used for IGa , IgG and Cd139 in the validation cohort\n- Figure 3A is not clear and C6 is not a clearly defined population.\n- Figure 3G HPV tumor appears to have a higher proportion of IGA in comparison to HPV + tumor contradicting the concluding message.\n- Figure 3 O and 3P please perform the tumor vs plasma Elisa from the same patient population to allow for definitive conclusion. If that is not possible include only the paired samples in the analysis.\nConclusion: need to be redirected to highlight the above limitations and not have misleading conclusions.\nOverall, the manuscript would benefit from description of the ethical consideration especially for tissue use of pts samples in the validation cohort.\nMinor comments:\nThe manuscript is missing important citations in support that are not included.", "reviewer_2": "In this manuscript, Tao et al. firstly established a comprehensive single-cell landscape of penile squamous cell carcinoma (PSCC), and found that there are a large number of immune cell infiltrations in PSCC. Compared to HPV-negative PSCC, HPV-positive PSCC patients exhibited significantly increased intratumor B cell, plasma cell, tertiary lymphoid structure (TLSs), IgA secretion and PIGR+ epithelial cells. Then, they used PSCC primary cells, xenograft tumor models and mouse genetic tools to demonstrate IgA-PIGR axis promotes anti-tumor activity. The higher expression of IgA-PIGR is associated with improved survival in PSCC. This study pioneered the establishment of a comprehensive cellular atlas of PSCC, shedding light on a novel anti-tumor response referred to as the B cell-IgA-PIGR axis within PSCC. However, the study still has numerous deficiencies in terms of insufficient evidence, conflicting viewpoints, and descriptive errors, which are as follows.\nMajor comments\n1. Based on the findings presented in Figures 2G, 2I, and 2J, there is a noteworthy and statistically significant increase in the number of tertiary lymphoid structures (TLS) within the tumor tissue. Some of these structures even approached a count of 1,000, while others exceeded 100 per square centimeter. These results are quite remarkable and warrant attention. Can you provide the original IHC images?\n2. Based on the findings presented in Figures 5T-W, it is evident that PIGR plays a role in mediating the inhibitory effects of IgA on tumor cell proliferation and metastasis, while also promoting apoptosis. However, the results in Figures 5C-E demonstrate that a similar effect can be achieved solely through the overexpression of PIGR. In this process, what role does IgA play?\n3. It is crucial to acknowledge and discuss the relevant literature that has not been cited. For example, PMID 31200747 demonstrates that B cell-deficient mice and B cell-depleted mice exhibit similar control of CRC compared to control mice. However, the findings in Figure 6B indicate that B cell deletion can actually inhibit MC38 tumor growth. This discrepancy should be addressed and analyzed to provide a comprehensive understanding of the observed results.\n4. The source of tumor IgA needs to be identified. Specifically, it is important to determine whether it is produced by tumor-infiltrating B cells. Notably, studies such as PMID 38447573, 37657842, and 37801341 demonstrate that there are very few B cells present in the subcutaneous tumor model, including the MC38 model. These references should be discussed to address the potential discrepancy and shed light on the cellular origins of tumor IgA.\nMinor comments\nFigure 1\n1. Why is there no sample “PN4”? please explain it clearly in the methods.\n2. According to the observations depicted in Figures 1K and 1L, it appears that there is no statistically significant variation in the extent of B cell and plasma cell infiltration between the tumor tissues and NAT tissues.\nFigure 2\n3. Does the infiltration of immune cells exhibit a higher prevalence in HPV-positive patients compared to HPV-negative patients?\n4. The textual description of Figures 2C, 2D, and 2E does not align with the visual representations depicted in the respective figures. It is necessary to rectify this inconsistency to ensure accuracy and coherence between the text and the figures.\nFigure 3\n5. In Figure 3A, there appears to be a discrepancy between the color representation of HPV- and HPV+ cells compared to the UMAP plot. The colors used in Figure 3A (Right) do not correspond accurately to the colors utilized in the UMAP plot.\n6. In Figure 3B, it is requested to present the statistical differences for each cell population, particularly for plasma cells, similar to the format demonstrated in Figure 1F. This would enable a comprehensive examination of the statistical variances among different cell populations.\n7. Figures 3G and 3H demonstrate that there were initially 4 HPV-negative patients and 4 HPV-positive patients who underwent BCR-Seq. However, according to the information provided in FS1A, it states that there were 2 HPV-negative patients (PT2, PT6) and 2 HPV-positive patients (PT16, PT19) involved in the BCR-Seq analysis. It is necessary to provide an explanation for this inconsistency to clarify the apparent discrepancy in the patient count.\nFigure 4\n8. In Figure 4A, the description stating \"6 paired NATs and PSCC tumors\" is inaccurate and fails to include the presence of \"PN4\".\n9. Regarding Figure 4K, it is important to clarify the rationale behind selecting proinflammatory cytokines and immunoglobulin as the underlying stimuli. The reasons for choosing these specific stimuli should be provided to elucidate the scientific basis for this experimental design.\n10. Figures 4N and 4O exhibit a discrepancy where the total sample number does not equal 90 as expected. It is necessary to explain the reason behind this inconsistency and provide a clear explanation regarding the actual sample count used in these figures.\nFigure 5\n11. In Figure S7B should say Re-sorting purity not Re-soringt purity.\n12. In Figures 5E and 5U should add image scale, including all IHC images. The image scale provides important information regarding the magnification or size of the images, ensuring accurate interpretation and comparison of the results.\n13. Is IgA colocalized with PIGR? Is this co-binding more prevalent in HPV-positive PSCC patients?\nFigure 6\n14. Figure 6N is inconsistent with the description in the text. The number of patients with double higher IgA and PIGR expression and double lower IgA and PIGR expression is inconsistent in Figures 6N and 6O.", "reviewer_4": "Using a cohort of a cohort of 90 patients with penile cancer, Tao et al have characterized the immune-environment of this disease, comparing HPV infection-dependent differences. The authors found that HPV+ penile cancers are heavily enriched in intra-tumoral B cells, plasma cells, and TLSs, with abundant IgA secretion (predictive of outcome) and expression of polymeric IgA receptor (PIGR) in tumor cells. Higher levels of IgA and PIGR expression are associated with better outcome, whereas the effects of IgA:PIGR interactions are characterized.\nThe manuscript recapitulates many of the concepts recently reported for human ovarian and endometrial cancer (and in NSCLC to a lesser extent), although those manuscripts are blatantly ignored (e.g., PMID: 33536615, 34949671, 37909039). These data are nevertheless relevant because: 1) Little is known about the ME of penile tumors; much less in a cohort of this size; 2) the results represent the first independent confirmation of mechanisms reported in other tumors by a single group (although ignored in the references); 3) there are significant differences between HPV+ and HPV- tumors for B cell infiltration and IgA production, but not for other immune determinants; 4) the authors made extensive use of single cell OMICS, with interesting differences across HPV-driven tumors; and 5) the unusual density of TLS in HPV+ tumors is interesting and could illustrate the immunogenicity of these tumors.\nThe strength of this manuscript is in the human data. The mouse data are questionable and need to be confirmed, or ignored.\nIn addition to adding the references previously ignored, the authors must address the following major comments:\n1-The distinction of B cells and plasma cells in Figure 1 seems to be based on expression of IGH1 and IGHA1 to define the latter. Independently of the possibility that plasma cells could produce other isotypes (e.g., IGHA2 or IGG4), isotype switched B cells can also express those genes. CD138 expression (as shown at the protein level in Figure 3) would be a definitory marker of plasma cells, while CD138- plasmoblasts could be defined as CD19+CD20-CD38+CD27-high cells. in any case, better markers are needed to establish these categories, because >10% of plasma cells in tumor beds (Fig1L) is unusually high.\n2-The structures shown in Figure 2 are indeed very strongly suggestive of TLS, but the presence of other components such as PNAd+ high endothelial venules and, ideally, CD21 follicular dendritic cells, would confirm this definition.\n3-The data shown in Fig5D-G are puzzling. It has been reported that PIGR provides an oncogenic signal in at least hepatocarcinoma by activating ribosomes through a poorly characterized mechanism (PMID: 33390805). If expression of PIGR is always deleterious for the tumor, even in the absence of IgA, why tumor cells do not evolve to lose PIGR?\n4-For Fig5J-R, previous studies reported RNA-seq changes elicited by IgA transcytosis through endometrial (PMID: 34949671) and ovarian (PMID: 33536615) cancer cells. How do these results compare to those analyses?\n5-The relevant data of this manuscript are the analyses of clinical samples, but: 1) why overexpressing Pigr in MC38, as opposed to CRISPR naturally expressed Pigr?; and 2) Do MC38 tumors really induce the production of IgA responses? This needs to be shown and it would be weird, because rapidly progressing mouse models induce limited isotype switched responses.\nMINOR COMMENTS\n-There is a typo in the title of Figure 5M", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Response to the reviewers’ comments on NCOMMS-24-16489\nOverview:\nWe thank the reviewers for their supportive comments and valuable suggestions that helped to\nimprove this manuscript significantly. Major alterations in text are indicated with red font in the\nrevised manuscript.\nWhile detailed responses to each point are provided in the following pages, we have summarized\nour responses to address the common concerns raised by the reviewers and the editor.\nSummary-1. PSCC Sample Selection and Numbers for Sequencing and Other Assays\nTo help all reviewers understand the structure of our study and why certain samples were chosen\nfor specific experiments, we include a detailed description of the sample selection criteria, the\nspecific assays or sequencing performed on each cohort, as below.\n(1) Single-cell RNA sequencing _ Discovery phase\nWe collected fresh tumor specimens (PT1-PT6) and paired normal adjacent tissue (NAT) (PN1-\nPN6) from 6 treatment-naïve PSCC patients (P1-P6; 4 HPV+ and 2 HPV-) for 3’scRNA-seq in\ndiscovery cohort 1. We also sorted CD45+ immune cells from additional 4 treatment-naïve PSCC\npatients (P7-P10; 2 HPV+ and 2 HPV) for 5’scRNA-seq and single-cell B-cell receptor sequencing\n(scBCR-seq) in discovery cohort 2. Note that the 3’scRNA-seq library construction for the PN4\nsample failed due to an excessive number of dead cells, as detailed in Revised Figure S1B.\n(2) Bulk BCR sequencing _ Discovery phase\nWe collected fresh tumor tissues from 8 PSCC patients (P11-P18; 4 HPV+ and 4 HPV-) for bulk\nBCR-sequencing to evaluate the types and numbers of immunoglobulin (Ig) isotype in PSCC.\nDetails are provided in Revised Figure S1B.\n(3) Flow cytometry analysis _ Validation phase\nWe harvested CD45+ immune cells from fresh tumor tissues of 20 PSCC patients (9 HPV+ and\nmicroenvironment (TME). This analysis was used to validate the relationship between tumor-\ninfiltrating B cells (TIL-Bs) and HPV status in PSCCs. Note that the selection of these samples\nwas based on the availability of fresh tumor tissues, without additional specific selection criteria.\nDetails are provided in Revised Figure S1C.\n(4) ELISA analysis _ Validation phase\nWe freshly collected plasma and paired tumor tissues from 17 PSCC patients (10 HPV+ and 7\nHPV-) to perform ELISA and measure IgA and IgG levels in both plasma and paired tumor tissues.\nNote that 13 patients provided paired plasma and tumor tissues, while the remaining 4 patients\nprovided only plasma. Details are provided in Revised Figure S1C.\n1\n(5) Immunohistochemistry (IHC) _ Validation phase\nWe obtained 90 formalin-fixed, paraffin-embedded (FFPE) PSCC tissue sections (38 HPV+ and\nstaining. The stains were for p16, CD3, CD20, CD138(SDC1), IgA, IgG, PIGR and IFNγ. Due to\ndepletion of some tissue sections, only 57 slices (29 HPV+ and 28 HPV-) were available for IFNγ\nstaining. Details are provided in Revised Figure S1C.\n(6) Multiplex immunofluorescence (mIF) staining _ Validation phase\nWe randomly selected FFPE PSCC tissue sections for multiplex immunofluorescence staining.\nTwo different marker sets were used: (CD3, CD20, CD38, Bcl6, DAPI) and (CD3, CD20, PNAD,\nCD21, DAPI) to observe tertiary lymphoid structure (TLS). Additionally, for immune\nmicroenvironment analysis, we used a third marker set (CD4, CD8, CD20, CD56, CD86, DAPI)\non the same 6 sections paired with samples used for 3’scRNA-seq. Details are provided in\nRevised Figure S1C.\nRevised Figure S1. Clinical information and sample selection of PSCC patients that\nsubjected to scRNA-seq and validation studies\n(B) Details of sequencing samples used for scRNA-seq (P1-P10) and bulk BCR-seq (P11-P18).\n(C) Overview of PSCC sample selection and numbers across validation assays, including flow\ncytometry (FCM), ELISA, IHC and mIF staining.\nTo address the reviewers' concerns regarding the use of p16INK4a (p16) immunohistochemistry\n(IHC) staining for HPV status definition, we provide a clear rationale, supported by relevant\nliterature, as below.\n(1) Rationale for using p16 IHC staining\nAlthough molecular methods like HPV PCR or in situ hybridization (ISH) are more specific and\nprovide direct evidence of HPV infection, p16 IHC staining offers several advantages for\ndetermining HPV status: (A) it is more cost-effective and readily available in pathology labs. (B)\np16 overexpression is strongly correlated with high-risk HPV infection, a relationship that is well-\nestablished in PSCC(PMID: 31243533)1 and HNSCC (PMID: 34099645; PMID: 31924475)2, 3. (C)\np16 IHC staining also aids in the classification and differential diagnosis of PSCC subtypes, such\nas distinguishing between low-grade verruciform and high-grade solid penile tumors, which may\nbe challenging to differentiate morphologically1. (D) Due to this, the International Society of\nUrological Pathology (ISUP) recommends p16 IHC staining for detecting HPV infection and for\naiding in the classification of penile tumors, highlighting its clinical utility (PMID: 29105175)4\n(These points have been incorporated into the revised Methods)\n(2) Criteria for calling p16 IHC positive\nA commonly accepted criterion for calling p16 IHC positive is when at least 70% of the tumor cells\nshow strong, diffuse nuclear, and cytoplasmic staining. This threshold is widely used in the context\nof HPV-related cancers, such as HNSCC5 and PSCC6. We have added this description in the\nrevised Methods\n(3) Validation of HPV status in scRNA-seq samples\nFor the 10 PSCC patients subjected to scRNA-seq, we performed both p16 IHC staining and\nHPV16/18 DNA in situ hybridization to confirm HPV status. The strong association between p16\noverexpression and high-risk HPV infection ensures that the classification of HPV+ and HPV-\ntumor tissues is rigorous and justified (Revised Figure S1A, S1D-S1E)\nsubjected to scRNA-seq and validation studies\n(A) The clinical details of PSCC patients are summarized in the table, including patient age, tumor\nsize, tumor stage, HPV16/18 DNA status (detected by in-situ hybridization), and p16 expression\n(detected by immunohistochemistry). P1-P6 represent 6 treatment-naïve PSCC patients from\nwhom we collected both surgical tumor specimens and paired normal adjacent tissue (NAT) for\nscRNA-seq (discovery cohort 1). P7-P10 represents 4 treatment-naïve PSCC patients where we\nsorted CD45+ immune cells from tumor specimens and performed scRNA-seq and paired single-\ncell B-cell receptor sequencing (scBCR-seq) (discovery cohort 2). All samples were derived from\nprimary tumor sites, with no any evidence of lymph node (LN) or visceral metastases. (D-E)\nRepresentative images showing detection of HPV status: (D) HPV16/18 DNA in situ hybridization\nand (E) p16 immunohistochemistry staining. The red triangles mark positive signals for HPV16/18\nDNA. P2 is representative of an HPV- patient, while P5 is representative of an HPV+ patient.\nThe reviewers raised an insightful question regarding to use precise markers to differentiate\nbetween various B cell subsets, particularly in distinguishing plasma cells from isotype-switched\nB cells. In response to the reviewers’ suggestion, we have reanalyzed the tumor-infiltrating B cells\n(TIL-Bs) using a more refined set of markers.\nBased on this reanalysis, we have categorized TIL-Bs into six clusters:\nnaïve B cells (C1, MS4A1+IGHD+IGHM+)\nEGR1- or EGR1+ memory B cells (C2-C3, MS4A1+IGHD-IGHM-CD24+CD27+)\ngerminal center B cells (C4, MS4A1+MME+)\nplasmablast (C5, MS4A1loCD38+SDC1-MKI67+)\nplasma cells (C6, MS4A1-CD38+SDC1+)\nAligning with the reviewers’ suggestions, we included CD138 (SDC1) as a well-established\nmarker for plasma cells, which typically lose expression of CD19 and CD20 as they mature. So,\nplasma cells are defined by the expression of CD138 (SDC1) along with high levels of CD38. This\nensures that we are accurately identifying mature plasma cells within the tumor microenvironment.\nBesides, to differentiate plasmablasts, we used a combination of CD20 (MS4A1), CD38, and\nMKI67 markers. Plasmablasts are defined as CD38+ but lack CD138 expression and are\nCD20loMKI67+. Moreover, we used additional markers like CD24, CD27, and EGR1 to distinguish\nbetween different memory B cell subsets, rather than relying solely on immunoglobulin heavy\nchain genes.\nWe have made these adjustments to our analysis and updated in Revised Figure 3A-3B\naccordingly to ensure that our classification of B cell subsets and plasma cells is accurate and\nreflects current best practices in the field. The updated data provide a more precise\ncharacterization of the B cell composition within the tumor microenvironment, addressing the\nconcern of unusually high plasma cell percentages.\nRevised Figure 3A-3B. scRNA-seq analysis of TIL-Bs in PSCC\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing), illustrating six clusters. The clusters are coloured by cell type (left), individual patient\n(middle) and HPV status (right). (B) Violin plots showing the normalized expression of selected\nmarker genes corresponding to each TIL-Bs cluster.\nTo address the reviewer’s concerns regarding the clonotyping strategy for defining B cell clones\nin our scBCR-seq analysis, we analyzed B cell clonality using the scRepertoire package (version\n1.3.5)7, which is widely used for this purpose8. We employed the combineBCR,\ncombineExpression, and clonalDiversity functions to determine clonal groups based on paired\nheavy and light chains.\nClonotypes were assigned only for cells where high-quality paired heavy and light chains were\nsequenced. This is essential for accurately defining B cell clones and avoiding misclassification.\nTo ensure a biologically meaningful and stringent clonotyping strategy, we utilized the CTstrict\nfunction within scRepertoire, which defines clones based on the following criteria: (1) V(D)J gene\nusage. Clones were defined by identifying B cells that share the same V gene and J gene usage\nin their immunoglobulin sequences. This method ensures that clonally related B cells have\nundergone the same recombination events during B cell development. (2) CDR3 sequence\nidentity. The CDR3 (complementarity-determining region 3), the most variable and functionally\nsignificant region of the B cell receptor, was compared across cells. To define clonally related\ncells, we used a >85% normalized hamming distance for the nucleotide sequence of the CDR3\nregion. This threshold ensures that minor variations due to sequencing errors or somatic\nhypermutation are accounted for, while maintaining high stringency for identifying clonal\nrelationships.\nBy employing this combination of V(D)J gene usage and CDR3 sequence identity, we ensured\nthat our clonotyping strategy was rigorous and aligned with best practices in the field. This method\ngoes beyond simply counting unique VH sequences, as it also includes key determinants such\nas J gene usage and CDR3 sequence similarity, which are critical for defining B cell clones with\nhigh confidence. Our methodology is consistent with recent studies that have employed similar\nstrategies to assess B cell clonality and diversity using scBCR-seq. We have updated the\nMethods section to more clearly articulate this approach.\nI would like to salute the authors efforts in conducting the study titled “A B cells-IgA-epithelial axis\nenhances antitumor response and improves outcome in HPV-associated penile squamous cell\ncarcinoma”. I have major clinical and methodological concerns for this study as it would benefit\nfrom answering the below major comments:\nWe appreciate the reviewer for the insightful questions and valuable suggestions that help to\nimprove this manuscript significantly. We have carefully reviewed each specific comment and\nprovided detailed responses to address the reviewer’s concerns.\nAbstract:\nQ1: Please clarify: “we collected a cohort of 90 PSCC patients and performed single-cell RNA\nanalyses on intratumoral cell populations from 10 treatment-naive primary PSCC cases.” Why\nonly 10 samples from the 90 PSCC pts collected were used any selection bias to clarify?\nResponse: We thank the reviewer for this question. The selection of 10 PSCC cases for scRNA-\nseq was based on the availability of fresh, high-quality tumor tissue necessary for this type of\nanalysis. After performing scRNA-seq, we extended our study to a broader cohort by using\nformalin-fixed, paraffin-embedded (FFPE) tissue slices from 90 previously collected PSCC cases.\nThese FFPE samples were used for immunohistochemistry (IHC) staining to validate the findings\nfrom the scRNA-seq analysis. For detailed sample information, please refer to Summary-1 in the\nresponse letter and Revised Figure S1B.\nQ2: Please avoid the use of words such as “massive” “remarkably enriched” without values to\nsupport.\nResponse: We thank the reviewer for this suggestion and apologize for any confusion caused by\nthe use of qualitative terms such as “massive” and “remarkably enriched” without accompanying\nquantitative data. We have carefully revised the manuscript to avoid the use of such words.\nIntroduction: would benefit from improved referencing and structing such as:\nQ3: Avoiding to use concluding statements in the introduction before the methodology and results\npresentation “We revealed that the intra-tumoral B cells (TIL-Bs), plasma cells, tertiary lymphoid\nstructures (TLSs) and IgA production were significantly enriched in HPV+ PSCC. Interestingly, a\nsubset of HPV+ PSCC malignant epithelial cells expressed polymeric IgA receptor (PIGR) to\nmediate IgA transcytosis, and elicited transcriptional changes that contributes to hindering\nmalignant progression. PIGR-overexpression promotes tumor cell apoptosis and enhances anti-\ntumor response in PSCC primary cells and in xenograft tumor models, while the phenotype was\nabolished in IgA deficient mice. In PSCC and some other tumor types that derived from mucosa,\nthe higher IgA and PIGR expression in tumor tissues correlated with improved survival outcomes.\nCollectively, our work firstly established the comprehensive picture of the PSCC ecosystem, and\ndemonstrated that the B cell-IgA-PIGR axis contributes to the clinical advantage in HPV+ PSCC.”\nResponse: We thank the reviewer for this valuable suggestion. We have revised the introduction\nto focus on the background and rationale, gaps in knowledge, objective of this study, while the\ndetailed findings were placed in the results and discussion sections.\nQ4: In the abstract it is stated that they performed scRNA seq on 10 samples in the methodology\ndiscovery cohort used is 6 samples only, please correct and clarify\nResponse: We thank the reviewer for this question and apologize for any confusion caused by\nthe description of the patient number. In our study, we performed scRNA-seq on a total of 10\nPSCC patients, but this involved two distinct discovery cohorts with slightly different sample\nnumbers and purposes:\n(1) Discovery Cohort 1 (3’ scRNA-seq): This cohort included 6 patients (P1-P6) from whom we\ncollected fresh tumor tissues and paired normal adjacent tissues (NATs). We performed 3’\nscRNA-seq to analyze both the tumor cells and the immune cells within these tissues. Specifically,\nwe successfully analyzed 5 NAT samples and 6 tumor samples. One NAT sample (PN4) could\nnot be included in the final analysis due to excessive dead cells, which resulted in poor-quality\ndata. This accounts for a total of 11 samples in this cohort.\n(2) Discovery Cohort 2 (5’ scRNA-seq and scBCR-seq): This cohort consisted of 4 additional\nPSCC patients (P7-P10) from whom we collected fresh tumor tissues and sorted CD45+ immune\ncells. For this cohort, we performed 5’ scRNA-seq and single-cell B-cell receptor sequencing\n(scBCR-seq) to focus specifically on the immune cells present in the tumor microenvironment.\nThis accounts for an additional 4 samples.\nWe have revised the manuscript and summarized the sample information in the response\nSummary-1 and Revised Figure S2B to ensure consistency and clarity regarding the number of\nsamples used and the specific methodologies applied to each cohort.\nRevised Figure S2B. Sample information for scRNA-seq data\n(B) Summary table of scRNA-seq data, detailing patient ID, tissue type (tumor or NAT), 10X library\nID, cell types, and the number of cells captured for each patient sample.\nfurthermore the second cohort is from the same institution? Need to clarify.\nResponse: We thank the reviewer for this question.\nAll PSCC cases used in our study, including both the discovery and validation cohorts, were\ncollected from the First Affiliated Hospital of the University of Science and Technology of China.\nWe confirm that both cohorts are from the same institution.\nIn our study, only fresh tumor tissues were suitable for scRNA-seq, which we applied to the\ndiscovery cohort of 10 PSCC patients. The validation cohort consisted of 90 PSCC cases with\nformalin-fixed, paraffin-embedded (FFPE) tissue slices. These samples were collected over the\npast 10 years and represent a larger clinical cohort for validating the findings from our discovery\ncohort at the protein and histological levels using immunohistochemistry (IHC). Due to the fixation\nprocess and long-term storage, FFPE samples are not suitable for scRNA-seq, which is why\nscRNA-seq was not applied to the validation cohort.\nWe have updated the manuscript to provide a clearer explanation of the study design.\nQ6: Ned to clarify the methodology for calling p16 IHC positive vs negative this is not clearly\ndetailed in the methodology section and is the main message of this manuscript. Multiple\nquantitative and qualitative methodology for p16 IHC calling can be used and this is a major\nimportant limitation for the validation cohort as no HPV PCR nor ISH was performed.\nResponse: We thank the reviewer for this question.\n(1) The methodology for calling p16INK4a (p16) IHC positive or negative: The determination of\np16 IHC positivity was based on established criteria commonly used in clinical practice and\nresearch for HPV-related cancers, including penile squamous cell carcinoma (PSCC) and\noropharyngeal squamous cell carcinoma (OPSCC). In our study, p16 positivity was defined as\nstrong and diffuse nuclear and cytoplasmic staining in at least 70% of the tumor cells. This\nthreshold is widely recognized and has been validated in multiple studies to correlate well with\nHPV-associated oncogenesis5, 6.\n(2) Rationale for using p16 IHC in the validation cohort: We acknowledge that molecular\nmethods like HPV PCR or in situ hybridization (ISH) offer more specific detection of HPV DNA or\nRNA, providing direct evidence of HPV infection. However, we used p16 IHC staining in the\nvalidation cohort for the following reasons: (A) p16 IHC is a cost-effective and widely available\nassay in pathology laboratories, making it a practical choice for routine clinical use and large-\nscale studies. (B) p16 overexpression is strongly associated with high-risk HPV infection. This\nrelationship is particularly well-established in PSCC1 and HNSCC2, 3, making p16 a reliable\nsurrogate marker. (C) p16 IHC staining is valuable in the classification and differential diagnosis\nof PSCC tumor subtypes. For example, it helps differentiate between low-grade verruciform and\nhigh-grade solid penile tumors, which can be challenging to distinguish morphologically1. (D) Due\nto this, the International Society of Urological Pathology (ISUP) recommends p16 IHC staining for\ndetecting HPV infection and aiding in the classification of penile tumors, underscoring its clinical\nutility4\nWe have revised the manuscript to clearly describe these points in revised Methods, including\nthe criteria used for p16 IHC positivity and the rationale for its use.\nmethodology reference.\nResponse: We thank the reviewer for this question and the opportunity to clarify the methodology\nused for staining and identifying tertiary lymphoid structures (TLS) in our study.\n(1) Methodology for staining and identifying TLS. In our study, we used CD3 and CD20\nimmunohistochemistry (IHC) to identify tertiary lymphoid structures (TLS) in our 90 FFPE PSCC\ntumor sections. CD3 and CD20 are well-established markers for T cells and B cells, which are\nthe primary components of TLS. The use of CD3/CD20 IHC staining is a widely accepted method\nfor the initial identification of TLS based on the presence and spatial organization of T and B cells.\nThis approach has been consistently employed in numerous studies to assess the immune\nlandscape within tumors9, 10, 11, 12.\nWhile CD3 and CD20 are sufficient for identifying TLS based on the organized clustering of T and\nB cells, we acknowledge that additional markers such as CD21 (for follicular dendritic cells) and\nPNAd (for high endothelial venules) are recommended to fully characterize the detailed\narchitecture and maturation of TLS. These additional markers are particularly useful when\nassessing the functional status of TLS or distinguishing between immature and mature TLS.\nHowever, given that the primary goal of our study was to compare the number of TLS between\nHPV- and HPV+ PSCC tumor sections, the use of CD3/CD20 IHC staining provided an appropriate\nand sufficient method for our specific research objectives.\n(2) Criteria for TLS identification and quantification. To ensure consistency and accuracy in\nour analysis, our criteria for TLS inclusion: (A) Structured organization: TLS were defined by a\nclear, organized structure, with B cells and T cells forming distinct zones. (B) Size threshold: only\nTLS consisting of more than 50 cells were included in the analysis to distinguish them from smaller,\nless organized immune cell clusters. (C) Normalization by tumor area: the number of TLS was\nnormalized by the area of the tumor sections to provide a consistent measure across samples.\nThis approach aligns with established methods for TLS identification and quantification in tissue\nsections, which often use similar size thresholds and normalization techniques to ensure robust\nand reliable data9, 13, 14.\nWe have revised the manuscript to clearly describe the methodology and criteria for TLS staining\nin revised Methods.\nResults:\nQ8: Figure 1 needs clarification on the methodology to highlight the number of cases and the\nlocation of the tumor collection penile primary alone or with LN and visceral metastatic sites.\nResponse: We thank the reviewer for this question.\nAll scRNA-seq tumor samples used in our study were collected exclusively from treatment-naïve\nPSCC patients (n=10). These samples did not include tissues from lymph node (LN) or visceral\nmetastatic sites. The selection criteria were strictly limited to primary penile tumor tissues to focus\non the intrinsic tumor microenvironment without the influence of metastatic spread.\nFor further clarity, the clinical information for the PSCC patients whose samples were used in the\nscRNA-seq analysis is detailed in the revised Figure S1A. This figure includes information on\nthe N and M stage of the tumors, where N0 and M0 indicates no lymph node involvement or\n10\ndistant metastasis. Therefore, all samples analyzed were from primary tumor sites without any\nevidence of LN or visceral metastases. We have updated the manuscript to explicitly state this\ninformation in the Revised Figure S1A legend.\nRevised Figure S1A: Clinical information of PSCC samples for scRNA-seq\n(A) The clinical characteristics of PSCC patients are presented in the table, detailing key\nparameters such as age, tumor size, tumor stage, HPV16/18 DNA in-situ hybridization (ISH)\nstatus, and p16 IHC results. Patients P1-P6 represent the 6 treatment-naïve PSCC cases from\nwhich surgical tumor specimens and paired normal adjacent tissues (NAT) were collected for\nscRNA-seq analysis (discovery cohort 1). Patients P7-P10 represent the 4 treatment-naïve PSCC\ncases in which CD45+ immune cells were isolated from tumor specimens for scRNA-seq and\npaired scBCR-seq (discovery cohort 2). All samples were collected from primary tumor sites, with\nno evidence of lymph node (N0) or distant metastasis (M0) as indicated by the N and M staging.\nQ9: The number of pts is not clear they state 6 treatment naïve PSCC in the results section.\nResponse: We thank the reviewer for this question. The total number of patients involved in the\nscRNA-seq analysis is 10. We have provided a detailed description of the sample selection and\ncohort composition in Summary-1 and Revised Figure S1. The clinical information is clearly and\nconsistently presented throughout the manuscript to avoid any further confusion.\nvalidation studies\n(B) Details of sequencing samples used for scRNA-seq (P1-P10) and bulk BCR-seq (P11-P18).\n(C) Overview of PSCC sample selection and numbers across validation assays, including flow\ncytometry (FCM), ELISA, IHC and mIF staining.\nQ10: Using HPV 16 and 18 DNa testing is not sufficient to detect all possible HPV subtypes, this\nhas potential major implications on the results as the paper is focused on the differences between\nHPV + and -.\nResponse: We thank the reviewer for highlighting the limitations associated with using only\nHPV16/18 DNA testing to detect HPV infection.\nFor the 10 PSCC patients subjected to scRNA-seq, we initially performed p16 IHC staining to\nassess HPV status. This approach aligns with the recommendations of the International Society\nof Urological Pathology (ISUP), which advocates for the use of p16 IHC staining to detect HPV\ninfection and assist in the classification of penile tumors due to its strong correlation with high-risk\nHPV infection (PMID: 31243533)1.\nTo further complement the p16 IHC results, we conducted HPV16/18 DNA ISH. We chose to\nfocus on HPV16 and HPV18 because they are the most prevalent high-risk HPV subtypes\nassociated with penile cancers (approximately 75% of HPV-positive penile cancers, Figure\nbelow)15. We recognize that relying solely on HPV16/18 DNA testing does not capture all possible\nhigh-risk HPV subtypes.\nHowever, the combination of p16 IHC staining and HPV16/18 DNA ISH enhances our ability to\naccurately classify HPV+ and HPV- tumor tissues. p16 IHC provides a broad marker for high-risk\nHPV-related oncogenesis, while HPV16/18 DNA ISH offers specific detection of the most\ncommon subtypes. Together, these methods provide a rigorous and justified approach to\ndetermining HPV status, especially when considering the prevalence of HPV16/18 in PSCC.\nrobustness of our findings (revised Methods and Revised Figure S1A, S1B-S1C)\nFigure: Pooled prevalence of HPV subtypes in penile cancer worldwide\n(A) The pooled prevalence of HPV6, HPV11, HPV16, HPV18, HPV31, HPV33, and HPV45 in\npenile cancer according to whether the cancer was positive for multiple HPV types or for a single\nHPV type, worldwide. Ranges are 95% CIs. HPV=human papillomavirus. (Lancet Oncol, 2019)\nQ11: Table 1 supplement has major issues: Missing HPV status per PCR only presented data on\np16 IHC, will need to include multiple methodology and to assess the discrepancy between p16\nIHC and HPV PCR or other HPV validated testing methods.\nResponse: We thank the reviewer for highlighting the need for clarity regarding the HPV status\ndata presented in Supplementary Table 1.\nIn Revised Figure S1A, we used both p16 IHC and HPV DNA ISH to cross-validate the HPV\nstatus of each sample. The consistent results from both methodologies (i.e., all p16-positive cases\nwere also positive for HPV16/18 DNA by ISH, and all p16-negative cases were negative for\nHPV16/18 DNA) provided strong validation for the HPV status of the 10 PSCC patients included\nin our scRNA-seq analysis.\nWhile p16 IHC and HPV DNA ISH are robust and widely used methods, we acknowledge that\nHPV PCR or other broad-spectrum HPV testing methods could detect additional high-risk HPV\nsubtypes beyond HPV16 and HPV18.\nWe appreciate the reviewer’s suggestion and have revised the supplementary materials and\nmanuscript to provide a clearer explanation of our methodologies and the rationale behind them.\nQ12: in table S1 the author has major errors that raise concern on the validity of the data:\nThey have 42 patients with positive tumors on Lymph node dissection so we would expect a minimum\nof 42 patient with N1 or N2 disease they only include 15 pts in the N stage. Also, they have 17\npts as stage 3 which requires N1 or N2 which would be an error with the 42 positive LND and also\nwith the 8 pts with N1 and the 9 pts with N2.\nResponse: We thank the reviewer for pointing out these discrepancies and for the opportunity to\nclarify and correct the data presented in Table S1. We apologize for the confusion caused by the\n13\nincorrect labeling in the original manuscript. The correct information is that 42 patients\nunderwent the lymph node dissection (LND) procedure, not that all had positive LND\nresults. This error led to the incorrect assumption that all 42 patients should have been\ncategorized as N1 or N2.\nAfter reviewing the data, we identified a statistical error in the categorization of patients by the\nAmerican Joint Committee on Cancer (AJCC) stages II and III. We have corrected the data to\nreflect that 54 patients are classified under AJCC stage II and 13 patients under stage AJCC III,\nrather than the previously stated 50 patients in stage II and 17 patients in stage III. We sincerely\napologize for these errors and any confusion they may have caused. To prevent similar issues in\nthe manuscript, we have implemented a more rigorous data verification process.\nQ13: Fig S1 a P1 T3N0M0 is placed as stage III and Pt19 t3n0M0 is stage II\nResponse: We thank the reviewer for identifying this error. We have corrected the staging\nclassification in Figure S1 to accurately reflect that P1 with T3N0M0 should be classified as AJCC\nstage II, not stage III. We apologize for this mistake. The manuscript has been updated\naccordingly to ensure consistency and accuracy.\ncohort and perform similar comprehensive testing. It is not feasible to make definitive conclusions\nsince the discovery cohort 1 only included 2 HPV negative pts.\nResponse: We sincerely appreciate the reviewer’s insightful suggestion regarding the structure\nof our discovery cohorts. We agree that an integrated analysis could provide a more\ncomprehensive understanding. However, several technical challenges prevented the successful\nintegration of our two scRNA-seq datasets:\n(1) Differences in library preparation and sequencing platforms: We used two different\nscRNA-seq platforms: 3' scRNA-seq for cohort 1 and 5' scRNA-seq for cohort 2. These methods\nare known to introduce variability, particularly in library preparation and sequencing depth (e.g.,\ndifferences in the number of reads per cell). Such variability can lead to inconsistencies in gene\nexpression profiles, making direct integration of these datasets prone to generating artifacts.\n(2) Cell composition variability between cohorts: Another critical difference was in the cell\nselection and dissociation protocols. In cohort 1, we captured the full spectrum of TME cells using\n3' scRNA-seq, while in cohort 2, we sorted only CD45+ immune cells post-dissociation, using 5'\nscRNA-seq. This led to an underrepresentation of myeloid populations in cohort 2. Myeloid cells,\nparticularly macrophages and dendritic cells, are sensitive to dissociation and sorting procedures,\noften resulting in cell death or loss during processing. Moreover, myeloid cells tend to adhere\ntightly to the extracellular matrix, complicating their recovery during dissociation. Consequently,\ncohort 2 displayed a skewed immune cell composition with an overrepresentation of lymphocytes,\nparticularly T and B cells, further complicating efforts to harmonize the two datasets.\nDespite these challenges, we still attempted to integrate the data. By isolating CD45+ immune\ncells from both cohorts, we clustered the immune cells into four main populations: TIL-Ts, TIL-Bs,\nmyeloid cells, and NK cells. Importantly, across both cohorts, we observed a consistent and\nsignificant increase in TIL-Bs in HPV+ PSCC compared to HPV- PSCC (Figure below). This\nconsistent finding, seen independently in both datasets, reinforces the robustness of our\nconclusion regarding the enrichment of TIL-Bs in HPV+ tumors.\nAlthough the integration was not ideal, this consistency across datasets supports the validity of\nour key conclusions. We have updated the manuscript to reflect this, emphasizing the robustness\nof the results despite the challenges faced with data integration.\n(A)UMAP plot of 119,351 single cells from the integrated scRNA-seq analysis of discovery cohort\n1 and 2. Cells are color-coded according to their classification into major cell types: immune cells,\nepithelial cells, endothelial cells, and cancer-associated fibroblasts (CAFs). (B) Immune cells\nwere identified based on the expression of the known immune-cell specific marker gene PTPRC\n(also known as CD45). Only cells expressing immune cell markers were bioinformatically isolated\nfor further analysis. (C) UMAP plot of 86,071 immune cells from both cohorts, categorized into\nfour distinct cell types: TIL-Ts, TIL-Bs, NK and myeloid cells. (D) The relative abundance of\nimmune cell types from each PSCC tumor, with patient labeled as P1-P10. The HPV status of\neach tumor is indicated above with ‘-’ for HPV-negative and ‘+’ for HPV-positive. (E) Bar graph\nshowing the relative proportions of immune cell types (TIL-Ts, TIL-Bs, NK cells, and myeloid cells)\nin HPV- (n=4) and HPV+ (n=6) PSCC tumors. Data are presented as mean ± s.e.m. The\nsignificance is indicated with *P<0.05.\nQ15: Fig 1G the used immune markers are not appropriate to evaluate the immune\nmicroenvironment, additional markers like Cd8, Cd4, Cd86, Cd263, cd168 and other immune\nmarkers of interest would have been more representative of the author’s message.\nResponse: We thank the reviewer for this comment regarding the evaluation of the immune\nmicroenvironment.\nincluding immune cells, malignant/non-malignant epithelial cells, endothelial cells, and cancer-\nassociated fibroblasts (CAFs) in normal adjacent tissues (NATs) and tumor tissues.\nIn response to the reviewer’s recommendation, we have expanded our analysis in the revised\nmanuscript. We performed additional staining using multiplex immunofluorescence (mIF) to\ninclude CD4, CD8, CD20, CD86, and CD56 markers, which better represent the diverse immune\ncell populations within the tumor microenvironment. These markers allow us to more precisely\nevaluate key immune cell subsets, such as:\n• CD4+ T cells (helper T cells)\n• CD8+ T cells (cytotoxic T cells)\n• CD86+ cells (indicative of antigen-presenting cells like dendritic cells and macrophages)\n• CD20+ B cells\n• CD56+ natural killer (NK) cells\nThe results, now presented in Revised Figure 1L-1M, demonstrate that immune cells expressing\nCD4, CD8, CD20, CD86, and CD56 are significantly increased in tumor tissues compared to\nNATs. This expanded analysis provides a more comprehensive understanding of the immune\nlandscape in PSCC, aligning with the reviewer’s suggestion and enhancing the relevance of our\nfindings regarding the immune microenvironment.\nRevised Figure 1L-1M. mIF staining of the tumor microenvironment in PSCC\n(L) Representative multiplex immunofluorescence images illustrating five major immune\npopulations in PSCC tumor and NATs. Shown are CD4+ T cells, CD8+ T cells, CD86+ antigen-\npresenting cells (APC, including macrophages and dendritic cell), CD20+ B cells, CD56+ natural\nkiller (NK) cells, and DAPI+ nuclei. Scale bar: 200μm. (M) Quantification of CD4+, CD8+, CD86+,\nCD20+ and CD56+ cells per tumor area in paired tumor and NATs regions from PSCC (n=6). Cell\nnumbers are presented as mean cell density per unit area.\ncohort 1, including 6,029 cells from HPV- tumor and 6,285 cells from HPV+ tumor (Fig. 2A)” Need\nto clarify the methodology for this sentence to select the immune cells as this cohort only has 2\nhpv negative vs 4 HPV + tumor samples, this raises concerns of selection bias and quality of data\nanalysis performed thereafter in that section of the manuscript.\nResponse: We thank the reviewer for raising this question regarding the methodology used for\nselecting immune cells and potential concerns about selection bias.\n(1) Clarification of the methodology for immune cell selection: To ensure the accurate\nselection of immune cells from the scRNA-seq datasets, we employed a rigorous bioinformatic\npipeline. (A) Quality control (QC). All scRNA-seq data underwent stringent QC checks to filter out\nlow-quality cells and doublets. (B) Normalization and clustering. After QC, the raw count data\nwere normalized, followed by dimensionality reduction. Cells were then clustered using\nunsupervised clustering algorithms to identify distinct cell populations without prior assumptions.\n(C) Identification of immune cells. Immune cells are identified based on the expression of the\nknown immune-cell specific marker gene PTPRC (also known as CD45). Only cells expressing\nimmune cell markers were bioinformatically isolated for further analysis. To ensure accuracy, we\nmanually validated the expression of additional immune-specific markers using heatmaps, feature\nplots, and violin plots to confirm that the clusters corresponded to expected immune cell types.\n(2) Addressing concerns of selection bias and sample representation: We understand the\nconcern regarding potential selection. Here are the steps we took to mitigate this. (A)\nStandardized criteria. The selection of immune cells was based on standardized criteria using\nwell-established immune cell markers (e.g. PTPRC), ensuring that the process was objective and\nindependent of the sample type (HPV-negative vs. HPV-positive). (B)Unbiased clustering. The\nclustering algorithms used for identifying cell populations were unsupervised and did not take\nsample identity into account, reducing the risk of bias in identifying immune cell populations. (C)\nData integration and batch effect correction. To address potential batch effects and inter-sample\nvariability, we integrated the datasets from the 2 HPV-negative and 4 HPV-positive tumors using\nfunction FindIntegrationAnchors and IntegrateData by Seurat’s integration workflow. These\nintegration techniques help to align cells from different samples into a common reference space,\nensuring that immune cell selection is not disproportionately influenced by the sample number or\ndifferences between individual tumors.\nIn summary, our approach to selecting immune cells from scRNA-seq data was systematic and\nunbiased, involving multiple layers of QC, objective clustering, and marker-based identification.\nThis methodology minimizes concerns about selection bias and ensures data quality. We have\nincluded a more detailed description of this process in the revised Methods.\nQ17: From the discovery cohort 2 of 90 pts the authors need to highlight the reason to select only\nResponse: We thank the reviewer for their question regarding the selection of PSCC samples\nfor flow cytometry validation in original Figure 2D-E. We have summarized all clinical samples\nused in our study in Summary-1 and Revised Figure S1C.\nDue to the low incidence of PSCC, the availability of fresh samples for research purpose in quite\nlimited. We do not have a specific selection criterion for these samples beyond their availability;\n18\nthe numbers of samples obtained depend on the incidence rate and the capacity to collect fresh\ntissues on a monthly basis.\nFor the validation of our scRNA-seq data, we utilized formalin-fixed, paraffin-embedded (FFPE)\nPSCC tissue sections. Over the past 10 years, we have collected 90 FFPE PSCC tissue sections\n(38 HPV+ and 52 HPV-) from the USTC-affiliated hospital, which were available for\nimmunohistochemistry (IHC) staining. This larger cohort of FFPE samples was used to broadly\nvalidate the presence and characteristics of tumor-infiltrating immune cells as observed in our\nscRNA-seq data.\nHowever, for flow cytometry analysis, we were constrained by the availability of fresh PSCC tumor\ntissues. We managed to obtain fresh tumor tissues from 20 PSCC patients (9 HPV+ and 11 HPV-),\nwhich were the only samples available to us at the time of the flow cytometry experiment. These\ntumor samples were dissociated into single-cell suspensions for flow cytometry to analyze\nimmune cell composition within the tumor microenvironment (TME). The flow cytometry analysis\nspecifically aimed to validate the relationship between tumor-infiltrating B cells (TIL-Bs) and HPV\nstatus in PSCC as suggested by our scRNA-seq findings.\nQ18: Clarify selection of the N=37 for fig 2G\nResponse: We thank the reviewer for this question. Out of the 90 FFPE PSCC tissue slides\navailable for our study, only 37 slides (N=37) contained both tumor and adjacent non-tumor\nregions within the same tissue section. These paired regions are essential for analyzing the\npresence and number of tertiary lymphoid structures (TLSs) in both the intra-tumor and adjacent\nnon-tumor areas, allowing for a direct comparison within the same tissue environment.\nWe have updated the manuscript to clearly state the selection process for these tissue slides to\nensure transparency and clarity.\nQ19: Fig 2H can be moved to supp\nResponse: We thank the reviewer for this suggestion. To enhance the readability and flow of the\nmanuscript, we have moved Figure 2H to the supplementary materials.\nQ20: Fig 2 I and J have errors and conflicting results with text please clarify.\nResponse: We thank the reviewer for highlighting the inconsistencies between the figures and\nthe text.\n(1) Clarification of Figures 2I and 2J: in the original manuscript, Figure 2I displayed the total\nnumber of TLSs, CD20+ B cells and CD3+ T cells within TLSs without normalization by tumor area,\nwhereas Figure 2J presented the normalized counts of these elements relative to the tumor area.\nThis led to confusion and conflicting interpretations in the text due to the lack of clarity about the\nnormalization process.\n(2) Corrections made: to address this issue and enhance clarity, we have removed Figure 2I\nfrom the manuscript, as the non-normalized data were less informative for comparative analysis.\nWe have retained the normalized data in what was Figure 2J, as it provides a more accurate\nrepresentation of the number of TLSs, CD20+ B cells and CD3+ T cells within TLSs relative to\ntumor area, which is essential for understanding their distribution in the tumor microenvironment.\nmanuscript.\nAdditionally, we have carefully reviewed and corrected any discrepancies between the text and\nFigure 2K to ensure consistency and accuracy in the presentation of our results.\nRevised Figure 2H, S3B-3C. Enrichment of tertiary lymphoid structures (TLSs) in HPV+\nPSCC\n(2H, S3B-3C) Quantification of the numbers of TLSs (2H), CD3+ T cells (S3B) and CD20+ B cells\n(S3C) within TLSs per tumor area in HPV− (n=52) and HPV+ (n=38) PSCC samples. Data\ndemonstrate a significant enrichment of TLSs and CD20+ B cells in HPV+ PSCC tumors compared\nto HPV−, while CD3+ T cells show no significant difference between the two groups.\nQ21: Please explain the selection of only 4 HPV- and HPV+ pts for bulk BCR seq, none of the\ncohort has 4 HPV- and 4 HPV+ pts?\nResponse: We thank the reviewer for the question regarding the selection of 4 HPV+ and 4 HPV-\npatients for the bulk B-cell receptor (BCR) sequencing experiment. We have summarized all\nclinical samples used in our study in Summary-1 and revised Figure S1B.\nFor the bulk BCR sequencing experiment, we specifically required fresh PSCC samples, as\nformalin-fixed, paraffin-embedded (FFPE) samples are not suitable for this type of sequencing\ndue to the degradation and cross-linking of DNA and RNA that occurs during the fixation process.\nAdditionally, the samples used for flow cytometry analysis were not available for bulk BCR\nsequencing because they were already processed and dissociated for cell sorting and immune\ncell composition analysis.\nWe acknowledge that this is a relatively small cohort. However, the low incidence of PSCC and\nthe specific requirements for fresh, high-quality tissue suitable for bulk BCR sequencing\nconstrained our sample size. The selection of these 8 samples was based solely on the availability\nof fresh samples at the time of the study, without any specific bias or selective criteria.\nDespite the small cohort, the bulk BCR sequencing of these 8 samples provided important insights\ninto the clonal diversity and B-cell repertoire in HPV+ and HPV- PSCC. These findings enhance\nour understanding of the immune microenvironment in PSCC, particularly regarding the role of B\ncells in HPV-associated versus non-associated tumors.\nWe have clarified these points in the revised manuscript to provide a better understanding of the\nsample selection process and its impact on our study.\nResponse: We thank the reviewer for pointing out the need to reference the stratification levels\nfor CD138 expression used in our study.\n(1) Clarification on CD138 stratification levels: The stratification of CD138 expression in our\nstudy was based on a semi-quantitative scoring system that combines staining intensity with the\npercentage of positive cells. This method is commonly used in histopathological analysis to\nprovide a more nuanced assessment of protein expression.\n(2) Scoring methodology: CD138-stained images were analyzed digitally using ImageJ software,\nand CD138 protein expression was manually scored based on both the intensity of staining and\nthe fraction of positive cells. The scoring system used was as follows:\n(0) Negative: No immune cells positive for CD138.\n(1+) Low: 1–5% of immune cells positive for CD138.\n(2+) Moderate: 6–20% of immune cells positive for CD138 with moderate or strong intensity, or\n6–100% of immune cells with weak intensity.\n(3+) High: 21–75% of immune cells positive for CD138 with moderate or strong intensity.\n(4+) Very High: 76–100% of immune cells positive for CD138 with moderate or strong intensity.\n(3) Supporting references: This scoring system is consistent with those used in previous studies\nthat assessed CD138+ plasma cell density in tissue sections (PMID: 23370224; PMID:\n26763251)16, 17. Furthermore, a similar approach to scoring protein expression in PSCC sections\nhas been described in PMID: 3339153918.\nWe have included these references in the revised manuscript to support the use of this\nstratification method.\nA23: Please add references and explain the different scoring system used for IgA, IgG and Cd139\nin the validation cohort\nResponse: We thank the reviewer for this question. We would like to clarify the scoring system\nused for CD138, IgA, and IgG in our study.\nFor CD138, we considered both the staining intensity and the fraction of positive cells, as CD138\nexpression varies widely in both aspects and is an established marker for plasma cells. Staining\nintensity reflects the overall protein level in individual cells, while the fraction of positive cells\nprovides information about how widespread the expression is within the tissue.\nIn contrast, for IgA and IgG, the scoring was based solely on the fraction of positive cells. This\ndecision was made because, in our study, IgA and IgG are secreted proteins, and the fraction of\npositive cells is a more direct indicator of the presence of these immunoglobulins in the tissue.\nThe intensity of staining for secreted proteins may be less indicative of biological relevance\ncompared to membrane-bound or intracellular proteins like CD138, where intensity correlates\nmore directly with the protein’s functional level within cells.\nThe scoring systems we applied are consistent with those used in previous studies. For CD138,\nthe combination of staining intensity and positive cell fraction has been commonly used in studies\nof plasma cell tumors, including ovarian cancer (PMID: 26763251). Similar fraction-based scoring\nfor immunoglobulins has also been applied in studies analyzing immunoglobulin expression in\ntumor tissues (PMID:23370224). These methods ensure that both the extent and the strength of\n21\nmarker expression are considered. We have updated the manuscript to include these references\nand provide a detailed explanation of the scoring systems used for IgA, IgG, and CD138 in\nrevised Methods.\nA24: Figure 3A is not clear and C6 is not a clearly defined population.\nResponse: We thank the reviewer for this question. As described in Summary-3, we have\nreanalyzed the tumor-infiltrating B cells (TIL-Bs) using a more refined set of markers. Based on\nthis reanalysis, we have categorized TIL-Bs into six clusters, including naïve B cells (C1,\nMS4A1+IGHD+IGHM+), EGR1- or EGR1+ memory B cells (C2-C3, MS4A1+IGHD-IGHM-\nCD24+CD27+), germinal center B cells (C4, MS4A1+MME+), plasmablast (C5,\nMS4A1loCD38+SDC1-MKI67+) and plasma cells (C6, MS4A1- CD38+SDC1+). We have made\nthese adjustments to our analysis and updated Revised Figure 3A-3B accordingly to ensure that\nour classification of B cell subsets and plasma cells is accurate and reflects current best practices\nin the field.\nRevised Figure 3A-3B. scRNA-seq analysis of TIL-Bs in PSCC\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing) illustrating six clusters. The clusters are colored by cell type (left), individual patient\n(middle) and HPV status (right). (B) Violin plots showing the normalized expression of selected\nmarker genes corresponding to each TIL-Bs cluster.\nQ25: Figure 3G HPV tumor appears to have a higher proportion of IGA in comparison to HPV+\ntumor contradicting the concluding message.\nResponse: We appreciate the reviewer’s insightful question regarding the original Figure 3G.\nHowever, this does not contradict our conclusion that IgA levels and diversity are increased in\nHPV+ PSCC compared to HPV- PSCC.\nThe key issue is how relative proportions were calculated. The proportions of IgA and other\nisotypes were based on their frequency within the total pool of BCR sequences. While we\nobserved a significant increase in the total number and diversity of BCR sequences in HPV+\nPSCC, this increase in absolute counts did not necessarily lead to a higher relative proportion for\neach isotype.\nHPV- tumors, IgA levels and diversity are indeed higher in HPV+ PSCC. This reflects an overall\nincrease in B-cell activity in HPV+ tumors, supporting our conclusion that IgA plays a more\nprominent role in the immune response in HPV+ PSCC compared to HPV- PSCC. We have\nupdated it in Revised Figure S4D\nQ26: Figure 3O and 3P please perform the tumor vs plasma Elisa from the same patient\npopulation to allow for definitive conclusion. If that is not possible include only the paired samples\nin the analysis.\nResponse: We thank the reviewer for the helpful suggestion to include only paired plasma and\ntumor tissue samples in the ELISA analysis to draw a more definitive conclusion.\n(1) Clarification of sample usage: In our study, we freshly collected plasma and paired tumor\ntissues from a cohort of 17 PSCC patients (10 HPV+ and 7 HPV-). Out of these, 13 patients\nprovided both plasma and paired tumor tissue samples, while the remaining 4 patients provided\nonly plasma samples. Details are provided in rebuttal Summary-1 and Revised Figure 3L-3M.\nFor the ELISA analysis shown in original Figures 3O and 3P, we measured IgA and IgG levels\nusing the 13 paired plasma and tumor tissue samples. This approach allowed us to directly\ncompare antibody levels between the tumor microenvironment and systemic circulation within the\nsame patients. By focusing on paired samples, we minimized inter-patient variability and ensured\na more reliable comparison of antibody distribution between tumor tissue and plasma.\n(2) Exclusion of unpaired plasma samples: In line with the reviewer’s recommendation, we\nhave revised the analysis to include only the 13 paired samples. By excluding the 4 unpaired\nplasma samples, we ensured that the tumor-to-plasma comparisons were more robust and based\nsolely on paired data, allowing for a more definitive conclusion regarding antibody levels in the\nsame patients.\nRevised Figure 3L-M. Elevated IgA levels in tumor tissues and plasma of HPV+ PSCC\n(L-M) Enzyme-linked immunosorbent assay (ELISA) was used to measure IgA and IgG levels in\ntumor tissues (L) and plasma (M) from freshly harvested samples of HPV- (n=7) and HPV+ (n=6)\nPSCC patients. Only paired samples (n=13) were included in the analysis to ensure direct\ncomparison of antibody levels between the tumor microenvironment and systemic circulation\nwithin the same patient population, minimizing inter-patient variability. Data are presented as\nMean ± s.e.m. Significance is indicated as *P < 0.05.\nconclusions. Overall, the manuscript would benefit from description of the ethical consideration\nespecially for tissue use of pts samples in the validation cohort.\nResponse: We appreciate the reviewer’s comments regarding the need to highlight limitations in\nthe conclusions and provide further clarification on the ethical considerations for the use of patient\nsamples in our validation cohort.\nWe would like to confirm that all patient samples used in this study, including the formalin-fixed,\nparaffin-embedded (FFPE) tumor sections from 90 PSCC patients, were collected in full\ncompliance with ethical guidelines and regulations.\nSpecifically, the use of these samples was approved by The First Affiliated Hospital of University\nof Science and Technology of China Review Board (2016KT67). Written informed consent was\nobtained from all patients prior to sample collection, ensuring that they were fully aware of the\npurpose of the study and how their samples would be used for research purposes. Additionally,\nall patient data were anonymized to protect patient confidentiality and privacy. The study was\nconducted in accordance with the Declaration of Helsinki to ensure responsible handling and\nuse of human samples.\nWe have included this clarification in the revised Methods section to outline the ethical\nframework under which the study was conducted.\nMinor comments:\nThe manuscript is missing important citations in support that are not included.\nResponse: We thank the reviewer for this suggestion. We have carefully revised the manuscript\nand cited all related refences to support our conclusion.\nIn this manuscript, Tao et al. firstly established a comprehensive single-cell landscape of penile\nsquamous cell carcinoma (PSCC), and found that there are a large number of immune cell\ninfiltrations in PSCC. Compared to HPV-negative PSCC, HPV-positive PSCC patients exhibited\nsignificantly increased intratumor B cell, plasma cell, tertiary lymphoid structure (TLSs), IgA\nsecretion and PIGR+ epithelial cells. Then, they used PSCC primary cells, xenograft tumor\nmodels and mouse genetic tools to demonstrate IgA-PIGR axis promotes anti-tumor activity. The\nhigher expression of IgA-PIGR is associated with improved survival in PSCC. This study\npioneered the establishment of a comprehensive cellular atlas of PSCC, shedding light on a novel\nanti-tumor response referred to as the B cell-IgA-PIGR axis within PSCC. However, the study still\nhas numerous deficiencies in terms of insufficient evidence, conflicting viewpoints, and descriptive\nerrors, which are as follows.\nResponse: We thank the reviewer for the thorough assessment and valuable suggestions, which\nhave greatly contributed to improving our manuscript. In response to the reviewer’s concerns\nregarding insufficient evidence, conflicting viewpoints, and descriptive errors, we have carefully\naddressed each issue as follows:\nMajor comments\nQ1: Based on the findings presented in Figures 2G, 2I, and 2J, there is a noteworthy and\nstatistically significant increase in the number of tertiary lymphoid structures (TLS) within the\ntumor tissue. Some of these structures even approached a count of 1,000, while others exceeded\n100 per square centimeter. These results are quite remarkable and warrant attention. Can you\nprovide the original IHC images?\nResponse: We appreciate the reviewer’s insightful question regarding the substantial increase in\nthe number of tertiary lymphoid structures (TLS) within tumor tissue, as observed in original\nFigures 2G, 2I, and 2J.\n(1) Clarification of TLS assessment: In our analysis, TLSs were quantified by counting the\nnumber of TLSs in each immunohistochemistry (IHC) section and normalizing these values by\nthe tumor area (measured in square centimeters). We only included TLSs composed of more than\n50 immune cells to maintain accuracy and consistency in our findings. This methodology aligns\nwith established approaches for TLS quantification in the literature9, 13, 14.\n(2) Clarification of TLS counts: In the original Figure 2I, the total number of TLSs was presented\nwithout normalization by tumor area, which accounts for the higher counts approaching 1,000. As\nreviewer #1 suggested, these non-normalized data were less informative for comparative\npurposes and we have removed Figure 2I. In the original Figure 2J (now Figure 2H), we present\nthe TLS count normalized by tumor area. The highest observed TLS density was 288 per cm²\n(2.88 per mm²), with an average of 36 per cm² (0.36 per mm²), offering a more accurate\nrepresentation for comparative analysis.\n(3) Providing IHC Images: In response to the reviewer’s request for original IHC images, we\nhave selected two representative images of HPV+ PSCC tumor sections, each showing TLSs\nhighlighted with red circles in an area of 1 mm². These images clearly demonstrate the abundance\nof TLSs per square millimeter and corroborate the statistical findings presented in the manuscript.\nWe trust these clarifications and the addition of representative images will address the reviewer’s\nconcerns.\nhighlighted with red circles in an area of 1 mm². (I) In the original Figure 2I, the total number of\nTLSs was presented without normalization by tumor area, which accounts for the higher counts\napproaching 1,000. We have removed this figure as reviewer#1 suggested. (H) TLS count\nnormalized by tumor area. The highest observed TLS density was 288 per cm² (2.88 per mm²),\nwith an average of 36 per cm² (0.36 per mm²). Data are presented as Mean ± s.e.m. Significance\nis indicated as **P < 0.01.\nQ2: Based on the findings presented in Figures 5T-W, it is evident that PIGR plays a role in\nmediating the inhibitory effects of IgA on tumor cell proliferation and metastasis, while also\npromoting apoptosis. However, the results in Figures 5C-E demonstrate that a similar effect can\nbe achieved solely through the overexpression of PIGR. In this process, what role does IgA play?\nResponse: We thank the reviewer for this insightful question.\n(1) Role of PIGR overexpression alone: As shown in Figures 5C-E, overexpression of PIGR\ncan lead to significant anti-tumor effects, including inhibition of tumor cell proliferation and colony\nformation, and decrease of cell migration, even in the absence of IgA. This suggests that PIGR\noverexpression alone can activate downstream signaling pathways, likely due to receptor\nclustering or constitutive activation at the cell membrane. These effects may overlap with the\nsignaling pathways typically triggered by the IgA-PIGR interaction, explaining the similar biological\noutcomes observed. However, it is important to note that PIGR overexpression, while effective,\nmay only partially replicate the full scope of the signaling network activated by IgA binding.\n(2) Specific role of IgA in PIGR-overexpressing cells: Despite the significant anti-tumor effects\nobserved with PIGR overexpression, our results demonstrate that IgA enhances these effects\nwhen present. For example, in PIGR-overexpressing PSCC primary cells, IgA binding further\namplifies the anti-tumor response by promoting a more significant increase in apoptosis (Figure\nbelow). This indicates that IgA binding to PIGR modulates and amplifies signaling pathways in\nways that PIGR overexpression alone does not fully achieve. In particular, IgA may serve as a\n26\nkey ligand that fine-tunes or enhances the downstream effects of PIGR signaling, maximizing the\nanti-tumor response.\n(3) IgA effects in cells with endogenous PIGR expression: To further investigate the specific\nrole of IgA in mediating anti-tumor effects, we conducted experiments using the A549 cell line,\nwhich naturally expresses high levels of PIGR. In these cells, IgA treatment significantly inhibited\nproliferation, reduced colony formation, and decreased cell migration in transwell assays\ncompared to PBS and IgG-treated groups (Revised Figure S9B-S9D). This provides additional\nevidence that IgA, when interacting with endogenously expressed PIGR, exerts a clear anti-tumor\neffect. These results emphasize that IgA-PIGR interactions play a critical role in tumor\nsuppression, even when PIGR is not artificially overexpressed.\nIn summary, while PIGR overexpression can independently activate certain anti-tumor pathways,\nour results indicate that the interaction between IgA and PIGR is essential for maximizing these\neffects. IgA serves as a crucial ligand that enhances PIGR-mediated inhibition of tumor\nproliferation, colony formation and migration, and promotes apoptosis more effectively than PIGR\noverexpression alone. Thus, the IgA-PIGR axis plays a critical role in driving the full anti-tumor\nresponse.\nFigure: Enhanced apoptosis in PIGR-overexpressing PSCC cells upon IgA treatment\n(A) Representative FACS analysis showing cell apoptosis in PIGR-overexpressing PSCC primary\ncells treated with1×PBS, IgG and IgA. The flow cytometry plots display staining for early-\napoptosis (SP, single positive for Annexin V) and late-apoptosis (DP, double positive for Annexin\nV and PI) populations. Data are presented as Mean ± s.e.m., with statistical significance indicated\nas follows: **P < 0.01, ****P < 0.0001.\n(B-D) Effects of IgA on PIGR-expressing A549 cells. (B) CCK8 proliferation assay showing that\nIgA (0.5 μg/ml) significantly inhibits the proliferation of PIGR-high A549 cells compared to PBS\nand IgG-treated controls. (C) Colony formation assay demonstrating reduced colony formation in\nIgA-treated A549 cells. Quantification shows a significant decrease in the number of colonies\ncompared to controls. (D) Transwell migration assay illustrating decreased migration in IgA-\ntreated A549 cells. Representative images are shown, and the quantification indicates a\nsignificant reduction in cell migration. Scale bar: 100 μm. Data are presented as Mean ± s.e.m.\nSignificance is indicated as *P < 0.05, **P < 0.01, ***P < 0.001.\nQ3: It is crucial to acknowledge and discuss the relevant literature that has not been cited. For\nexample, PMID 31200747 demonstrates that B cell-deficient mice and B cell-depleted mice\nexhibit similar control of CRC compared to control mice. However, the findings in Figure 6B\nindicate that B cell deletion can actually inhibit MC38 tumor growth. This discrepancy should be\naddressed and analyzed to provide a comprehensive understanding of the observed results.\nResponse:\nWe thank the reviewer for highlighting this important point. We have reviewed the relevant\nliterature, including the study mentioned by the reviewer (PMID 31200747), and have cited the\nrelated references in our revised manuscript.\n(1) Differences in experimental timing and conditions: In the study cited by the reviewer\n(PMID 31200747), anti-CD20 antibodies were used to reduce tumor-infiltrating B cells in the\nMC38 colorectal cancer model. They initiated anti-CD20 antibody treatment on day 9, which was\nthree days after the start of anti-PD-1 therapy. Their results showed that B cell depletion under\nthese conditions did not affect tumor growth. In contrast, our study administered anti-CD20\nantibodies at the time of MC38 tumor cell implantation. We observed that B cell deletion at this\nearlier time point significantly inhibited tumor growth. This suggests that the timing of B cell\ndepletion may be a critical factor in its efficacy; early intervention might prevent the establishment\nor expansion of a pro-tumorigenic B cell population, whereas later depletion may not reverse\nestablished immune interactions within the tumor microenvironment.\n(2) Supporting evidence from other studies: Our findings are supported by another study\npublished in Cancer Cell (PMID: 24909985), where mice were treated with anti-CD20 monoclonal\nantibodies (mAb) either at the time of SCC tumor cell implantation and again at 2-week intervals\n(day0; day14), or after orthotopic tumors had reached ~50 mm³ (day12) and then day14 (Figure\nA below)19. This study found that administering anti-CD20 mAb concurrent with SCC implantation\nsignificantly slowed tumor growth to the end stage. However, when anti-CD20 mAb therapy was\ninitiated after tumors became apparent, the treatment did not significantly slow tumor growth. This\nfurther supports the hypothesis that early depletion of B cells may have a different impact on\ntumor progression compared to depletion at later stages.\n(3) Additional evidence from B cell-deficient models: Furthermore, a recent study published\nin Nature (PMID: 34732892) demonstrated that B cell-deficient mice (muMt-/-) showed significantly\ninhibited MC38 tumor growth compared to wild-type (WT) mice in placebo-treated control groups\n(Figure B below). This suggests that the absence of B cells can indeed inhibit MC38 tumor growth,\naligning with our observations that early B cell depletion via anti-CD20 treatment can suppress\ntumor growth.\nIn summary, the discrepancies between our findings and those of the study mentioned by the\nreviewer likely stem from differences in the timing of B cell depletion, the tumor stage at which\ntreatment is initiated, and potentially the tumor microenvironment's composition and dynamics at\nthose times. Early depletion of B cells might prevent their pro-tumorigenic functions or the\nestablishment of an immune-suppressive environment, whereas later depletion might not be as\neffective. These differences highlight the importance of considering the timing and context of\nimmune cell depletion strategies in cancer therapy. We will incorporate these considerations and\nreferences into the revised manuscript to provide a comprehensive understanding of the\nobserved results.\nin PDSC5.6-derived tumors in mice treated with anti-CD20 (aCD20) or control (aRW) mAb. Mice\ntreated with aCD20 from tumor inoculation (day 0, red) had significantly slower tumor growth\ncompared to those treated at day 12 (green) after tumors became visible. Statistical significance\nbetween groups is denoted by * (P < 0.05). Data shown as mean ± SEM, n = 10 per group. (B)\nMC38 tumor growth in WT and B cell-deficient mice. MC38 tumor growth in wild-type (WT)\nand B cell-deficient (muMt-/-) mice implanted with GABA or placebo pellets. Tumor growth was\nsignificantly inhibited in muMt-/- mice, particularly in placebo-treated groups. Data shown as mean\n± SEM, n = 6-8 per group.\nQ4: The source of tumor IgA needs to be identified. Specifically, it is important to determine\nwhether it is produced by tumor-infiltrating B cells. Notably, studies such as PMID 38447573,\n37657842, and 37801341 demonstrate that there are very few B cells present in the\nsubcutaneous tumor model, including the MC38 model. These references should be discussed\nto address the potential discrepancy and shed light on the cellular origins of tumor IgA.\nResponse: We thank the reviewer for this important question regarding the source of IgA in the\nMC38 tumor model.\n(1) IgA source in the MC38 model: We agree with the reviewer that studies such as PMID\n38447573, 37657842, and 37801341 highlight the limited presence of B cells in subcutaneous\ntumor models, including the MC38 model. This raises legitimate questions about the origin of IgA\ndetected in these tumors, suggesting that it is unlikely to be produced locally by tumor-infiltrating\nB cells or plasma cells.\nGiven the low infiltration of B cells and plasma cells in the MC38 tumor microenvironment, we\npropose the following potential sources for the IgA observed: (A) Plasma cells from draining lymph\nnodes: IgA may be produced by plasma cells that have differentiated from B cells in nearby\nlymphoid tissues, such as draining lymph nodes. These plasma cells can either migrate into the\ntumor microenvironment or secrete IgA that is then transported into the tumor, contributing to local\nIgA levels20, 21. (B) Mucosal-Associated Lymphoid Tissues (MALT): Another possible source is\nmucosal-associated lymphoid tissues (MALT), which are highly efficient at producing IgA. In\ncases of systemic immune responses, IgA generated in MALT could enter the tumor\nmicroenvironment through the circulatory or lymphatic system20, 21. (C) Systemic circulation and\nPIGR transcytosis: IgA could enter the tumor microenvironment from systemic circulation22. In this\nscenario, IgA might be transcytosed into the tumor through the Polymeric Immunoglobulin\nReceptor (PIGR), which is expressed by epithelial or tumor cells. This mechanism would explain\nthe presence of IgA in the tumor, even with minimal local B cell infiltration\n(2) Limitations of the MC38 model for local IgA production: We acknowledge that the MC38\nmodel has limitations for studying the direct local production of IgA, given the scarcity of B cells\nand plasma cells in the tumor microenvironment. However, our primary focus in this study was to\ninvestigate the downstream effects of IgA, particularly in relation to its interaction with PIGR\nexpressed on epithelial or tumor cells. The MC38 model still provides valuable insights for this\npurpose, especially considering its immunocompetent background and relevance to cancer\nimmunotherapy research.\n(3) Significance of the IgA-PIGR axis: Our findings underscore the importance of the IgA-PIGR\naxis in tumor suppression. Notably, in our study, the absence of IgA alone did not affect MC38\ntumor growth, suggesting that IgA may not be the primary driver of anti-tumor activity in this model.\nHowever, we observed enhanced anti-tumor effects with PIGR overexpression, and a significant\nreduction in these effects in the absence of IgA, indicating that an intact IgA-PIGR axis is crucial\nfor tumor suppression (Revised Figure 6H). These results highlight the significance of IgA-PIGR\ninteractions in modulating the tumor microenvironment.\nWe have revised the manuscript to incorporate these discussions, addressing the potential\nsources of IgA in the MC38 model and the broader context of our findings.\nRevised Figure 6H. Tumor growth curves in WT (IgA+/+) and IgA deficient (IgA-/-) mice inoculated\nsubcutaneously with 5´105 control- and PIGR- overexpressed MC38 cells. n=5. Data shown as\nmean ± SEM. Significance is indicated as *P < 0.05, **P < 0.01.\nMinor comments\nFigure 1\nQ5: Why is there no sample “PN4”? please explain it clearly in the methods.\na technical issue encountered during the scRNA-seq library construction. Specifically, the library\nconstruction for the adjacent normal tissue sample from patient P4 (PN4) failed due to an\nexcessive number of dead cells, which compromised the quality of the sample. As a result, we\nwere unable to include PN4 in the final analysis.\nWe have summarized the details of sample selection and the number of samples in Summary-1\nand Revised Figure. S1B. We will ensure this clarification is included in the methods section to\nprovide a comprehensive understanding of the sample availability and selection process.\nQ6: According to the observations depicted in Figures 1K and 1L, it appears that there is no\nstatistically significant variation in the extent of B cell and plasma cell infiltration between the tumor\ntissues and NAT tissues.\nResponse: We thank the reviewer for the question regarding the original Figures 1K and 1L,\nwhere the relative proportions of B cells and plasma cells in the scRNA-seq data show no\nstatistically significant differences between tumor tissues and adjacent normal tissues (NAT).\n(1) Limitations of scRNA-seq data: We acknowledge that our scRNA-seq analysis was based\non a relatively small cohort of PSCC samples. This limited sample size may affect the statistical\npower needed to detect significant differences in immune cell proportions between tumor and\nNAT tissues. Despite the lack of statistical significance in the relative proportions of B cells and\nplasma cells, we observed a clear trend of increased proportions of these cell types in tumor\ntissues, suggesting a biologically relevant pattern that may not have reached statistical\nsignificance due to sample size constraints.\n(2) Complementary validation through mIF: More importantly, while the scRNA-seq data did\nnot reveal statistically significant differences in relative proportions, our multiplex\nimmunofluorescence (mIF) analysis provided additional evidence. Specifically, the absolute\nnumber of B cells per tumor area was significantly increased in PSCC tumor tissues compared to\nNATs (Revised Figure 1L-1M). This finding supports our overall conclusion that PSCC tumors\nexhibit a higher infiltration of immune cells, including B cells, validating the trend observed in the\nscRNA-seq data.\nIn summary, we have revised the manuscript to present a more balanced interpretation of the\ndata. While we acknowledge the absence of statistical significance in the scRNA-seq analysis of\nrelative proportions, we emphasize the complementary results from mIF, which confirm a\nsignificant increase in B cells in tumor tissues. This additional evidence strengthens the\nconclusion of enhanced immune cell infiltration in PSCC tumors.\nFigure 2\nQ7: Does the infiltration of immune cells exhibit a higher prevalence in HPV-positive patients\ncompared to HPV-negative patients?\nResponse: We thank the reviewer for this important question regarding the prevalence of\nimmune cell infiltration in HPV-positive versus HPV-negative PSCC patients.\nTo assess this, we conducted flow cytometry (FCM) analysis on fresh tumor specimens from\nPSCC patients to examine the levels of tumor-infiltrating B cells (TIL-Bs) and tumor-infiltrating T\ncells (TIL-Ts). These immune cells were gated from CD45+ cells, representing the total immune\ncell population, as shown in the original Figures 2C-2D.\nFollowing the reviewer’s suggestion, we extended our analysis by evaluating the percentage of\ntotal CD45+ immune cells within the live cell population in both HPV- and HPV+ PSCC tumor\ntissues. Our results showed that there was no statistically significant difference in the overall\nproportion of CD45+ immune cells between HPV- and HPV+ tumor tissues (Revised Figure 2C-\nD). We have updated the manuscript to reflect these findings, providing a more nuanced\nunderstanding of immune cell infiltration in HPV- versus HPV+ PSCC tumors.\nRevised Figure 2C-2D. Flow cytometry analysis of TIL-Ts, TIL-Bs in HPV+ PSCC vs. HPV-\nPSCC.\n(C) Gating strategy used for flow cytometry (FACS) analysis to identify tumor-infiltrating B cells\n(TIL-Bs) and tumor-infiltrating T cells (TIL-Ts) in fresh tumor specimens from PSCC patients.\n(D) Histograms displaying the relative proportions of CD45+ immune cells, TIL-Bs, total TIL-Ts,\nCD4+ TIL-Ts, and CD8+ TIL-Ts in fresh tumor specimens of HPV− (n=11) and HPV+ (n=9) PSCC\npatients. Data are presented as mean ± s.e.m., showing no statistically significant difference in the\noverall proportion of CD45+ immune cells between HPV− and HPV+ tumors.\nQ8: The textual description of Figures 2C, 2D, and 2E does not align with the visual\nrepresentations depicted in the respective figures. It is necessary to rectify this inconsistency to\nensure accuracy and coherence between the text and the figures.\nResponse: We thank the reviewer for pointing out the inconsistency between the textual\ndescriptions and the visual representations in original Figures 2C, 2D, and 2E. We have carefully\nreviewed and corrected the text in the manuscript to ensure it accurately aligns with the data\ndepicted in these figures.\nQ9: In Figure 3A, there appears to be a discrepancy between the color representation of HPV-\nand HPV+ cells compared to the UMAP plot. The colors used in Figure 3A (Right) do not\ncorrespond accurately to the colors utilized in the UMAP plot.\nResponse: We thank the reviewer for highlighting the discrepancy in color representation\nbetween the UMAP plot and the dot plot in original Figure 3A. We have reviewed the figure and\ncorrected the colors used in Revised Figure 3A (Right) to ensure they accurately correspond to\nthe colors used in the UMAP plot.\nRevised Figure 3A. scRNA-seq analysis of TIL-Bs in PSCC\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing) illustrate six clusters. The clusters are colored by cell type (left), individual patient\n(middle) and HPV status (right).\nQ10. In Figure 3B, it is requested to present the statistical differences for each cell population,\nparticularly for plasma cells, similar to the format demonstrated in Figure 1F. This would enable\na comprehensive examination of the statistical variances among different cell populations.\nResponse: We thank the reviewer for the insightful suggestion to present statistical differences\nfor each cell population, particularly plasma cells, as demonstrated in original Figure 1F.\n(1) Reanalysis of TIL-Bs: In response to the reviewers’ request, we reanalyzed the TIL-B\npopulations using refined markers. This allowed us to classify the TIL-Bs into six distinct clusters:\nnaïve B cells (C1, MS4A1+IGHD+IGHM+), EGR1- or EGR1+ memory B cells (C2-C3,\nMS4A1+IGHD-IGHM-CD24+CD27+), germinal center B cells (C4, MS4A1+MME+), plasmablasts\n(C5, MS4A1loCD38+SDC1-MKI67+), and plasma cells (C6, MS4A1- CD38+SDC1+). These\nrefinements are now reflected in Revised Figures 3A-3C, ensuring precise classification aligned\nwith current best practices.\n(2) Statistical power and observed trends: Although our scRNA-seq analysis was limited in\ndetecting statistically significant differences between HPV+ and HPV- tumors due to the relatively\nsmall sample size and inherent variability in human samples, we did observe a clear trend toward\nan increased plasma cell population (C6) in HPV+ tumors (Revised Figures 3C). While this trend\ndid not reach statistical significance in the scRNA-seq data, it strongly suggested a potential\ndifference in plasma cell infiltration between HPV+ and HPV- tumors.\n33\n(3) Validation with IHC staining: To further validate the observed trends, we performed IHC\nstaining for CD138 to quantify plasma cell infiltration in both HPV+ and HPV- tumors. This analysis\nconfirmed a statistically significant increase in plasma cell infiltration in HPV+ tumors (Revised\nFigures 3F, 3I), consistent with the trends seen in our scRNA-seq data. By complementing the\nscRNA-seq data with IHC validation, we were able to provide stronger evidence supporting the\nincreased presence of plasma cells in HPV+ PSCC.\nRevised Figure 3. Plasma cells are increased in HPV+ PSCC tumor compared to HPV-\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing), illustrating six clusters. The clusters are colored by cell type (left). (C) Bar graph\ndisplaying the relative proportions of the TIL-Bs clusters, stratified by HPV status (HPV+ vs. HPV-)\nin PSCC tumors. (F) Bar graph representing the distribution of IHC scores for CD138 across HPV-\n(n=52) and HPV+ (n=38) PSCC tumor sections. IHC scores range from 0 to 4, with higher scores\nindicating greater CD138 expression. (I) Representative images and the mean IHC scores for\nCD138 staining in HPV- (n=52) and HPV+ (n=38) PSCC tumor samples. Scale: 100μm.\nQ11. Figures 3G and 3H demonstrate that there were initially 4 HPV-negative patients and 4 HPV-\npositive patients who underwent BCR-Seq. However, according to the information provided in\nFS1A, it states that there were 2 HPV-negative patients (PT2, PT6) and 2 HPV-positive patients\n(PT16, PT19) involved in the BCR-Seq analysis. It is necessary to provide an explanation for this\ninconsistency to clarify the apparent discrepancy in the patient count cell populations.\nResponse: We thank the reviewer for the question and apologize for any confusion regarding\nthe patient numbers.\nIn original Figure S1A, we performed 5' scRNA-seq coupled with scBCR-seq on CD45+ immune\ncells from a cohort of 4 treatment-naïve PSCC patients (now labeled as P7-P10), including 2 HPV-\n(P7, P8) and 2 HPV+ (P9, P10) patients. The scBCR-seq was designed to evaluate BCR clones\nat a single-cell resolution.\nIn the original Figures 3G and 3H, the data presented was from bulk BCR-seq analysis performed\non tumor tissues from 4 HPV- (now labeled as P11-P14) and 4 HPV+ (P15-P18) PSCC patients.\nThis bulk BCR-seq aimed to evaluate the overall diversity and composition of immunoglobulin (Ig)\nisotypes in PSCC tumors. These patients were distinct from those involved in the scBCR-seq\nanalysis.\nand another for scBCR-seq. We have updated the manuscript to explain the sample numbers\nused in Revised Figure S1A-C to ensure transparency and consistency.\nRevised Figure S1. Clinical information and sample number of PSCC patients that\nsubjected to scRNA-seq and validation studies\n(A) The clinical characteristics of PSCC patients are presented in the table. Patients P1-P6\nrepresent the 6 treatment-naïve PSCC cases from which surgical tumor specimens and paired\nnormal adjacent tissues (NAT) were collected for scRNA-seq analysis (discovery cohort 1).\nPatients P7-P10 represent the 4 treatment-naïve PSCC cases in which CD45+ immune cells were\nisolated from tumor specimens for scRNA-seq and paired scBCR-seq (discovery cohort 2). All\nsamples were collected from primary tumor sites, with no evidence of lymph node (N0) or distant\nmetastasis (M0) as indicated by the N and M staging.\n(B) Details of sequencing samples used for scRNA-seq (P1-P10) and bulk BCR-seq (P11-P18).\n(C) Overview of PSCC sample selection and numbers across validation assays, including flow\ncytometry (FCM), ELISA, IHC and mIF staining.\nQ12: In Figure 4A, the description stating \"6 paired NATs and PSCC tumors\" is inaccurate and\nfails to include the presence of \"PN4\".\nResponse: We thank the reviewer for pointing out this inaccuracy in the description of original\nFigure 4A. We apologize for any confusion caused by this omission. We have now updated the\nfigure legend to accurately reflect the actual number of samples analyzed. The revised Figure\n4A legend now reads: \"Heatmap of the large-scale CNVs profile for individual epithelial cells\n(rows) from 5 paired NATs and PSCC tumors (excluding PN4 due to failed library construction),\nusing the same scRNA-seq dataset as in Fig. 1A.\"\nQ13: Regarding Figure 4K, it is important to clarify the rationale behind selecting proinflammatory\ncytokines and immunoglobulin as the underlying stimuli. The reasons for choosing these specific\nstimuli should be provided to elucidate the scientific basis for this experimental design.\nResponse: We thank the reviewer for this insightful question.\nThe selection of proinflammatory cytokines and immunoglobulin as stimuli in Figure 4K was\ngrounded in both prior research and our own experimental findings. Previous studies have\ndemonstrated that certain proinflammatory cytokines, including interferon-gamma (IFN-γ), tumor\nnecrosis factor (TNF), and interleukin-1 (IL-1), are key regulators of Polymeric Immunoglobulin\nReceptor (PIGR) expression23 (PMID: 16048543). These cytokines play critical roles in immune\nactivation and inflammation, processes that are particularly relevant to the tumor\nmicroenvironment.\nIn our study, HPV+ PSCC-enriched clusters C3 and C4 displayed high PIGR expression and were\ncharacterized by genes related to immune activation and inflammatory responses (Revised\nFigure 4E). We also observed that both plasma cells and PIGR+ epithelial cells were significantly\nmore prevalent in HPV+ tumors. A positive correlation between PIGR and IgA levels in IHC\nanalysis (Revised Figure 4S) suggested a potential link between IgA and PIGR expression in\nthe tumor microenvironment.\nBased on these findings, we selected proinflammatory cytokines (IFN-γ, TNF, IL-1) and\nimmunoglobulin (IgA) as stimuli to investigate their role in upregulating PIGR. This approach\naimed to explore the mechanisms behind PIGR upregulation in HPV+ tumors and its potential\nrelationship with immune activity. We have updated the manuscript to clarify this rationale and\nbetween\nand better explain the choice of these stimuli.\nIFNγ\n(E) Gene Ontology (GO) analysis highlighting the upregulated differentially expressed genes\n(DEGs) in HPV+ PSCC-enriched clusters 3 and 4, showing significant enrichment in immune\nactivation and inflammatory processes, including innate immune responses and type II interferon\nsignaling. (S) Pearson correlation analysis of IHC scores demonstrating a positive correlation\nbetween PIGR and IgA expression in PSCC tumor sections,\nQ14: Figures 4N and 4O exhibit a discrepancy where the total sample number does not equal 90\nas expected. It is necessary to explain the reason behind this inconsistency and provide a clear\nexplanation regarding the actual sample count used in these figures.\nResponse: We thank the reviewer for pointing out the discrepancy in sample numbers in original\nFigures 4N and 4O.\nIn our validation cohort, we initially included 90 PSCC cases for various analyses, such as\nstaining for p16, CD3, CD20, CD138, IgA, IgG and PIGR. However, for IFNγ staining, only 57\nformalin-fixed, paraffin-embedded (FFPE) tumor sections (29 HPV+ and 28 HPV-) were available.\nThis reduction occurred due to the depletion of some tissue sections during previous staining\nprocedures.\nGiven that multiple analyses were performed on the same cohort, some tissue sections were\nexhausted, limiting the available samples for IFNγ staining. To ensure consistency and avoid\nintroducing variability, we did not add new tissue sections from other patients. Therefore, we\nproceeded with the 57 available samples.\nWe have clarified this in the Revised Figure 4 legends and provided detailed sample information\nin Summary-1 and Revised Figure S1C to ensure transparency regarding the number of\nsamples used in each analysis.\nFigure 5\nQ15: In Figure S7B should say Re-sorting purity not Re-soringt purity.\nResponse: We thank the reviewer for catching this typo. We have corrected the spelling from\n\"Re-soringt purity\" to \"Re-sorting purity\" in the revised manuscript.\nQ16：In Figures 5E and 5U should add image scale, including all IHC images. The image scale\nprovides important information regarding the magnification or size of the images, ensuring\naccurate interpretation and comparison of the results.\nResponse: We thank the reviewer for this suggestion. We have added images scale for all IHC\nimages, which helps improve the clarity and accuracy of our presentation.\nQ17: Is IgA colocalized with PIGR? Is this co-binding more prevalent in HPV-positive PSCC\npatients?\nResponse: We thank the reviewer for this insightful question.\n(mIF) staining to assess the colocalization of IgA and PIGR. Our analysis revealed that IgA and\nPIGR do colocalize within the tumor tissue, although not all IgA is co-localized with PIGR.\nMoreover, we observed a significant increase in the overall number of co-localized IgA and PIGR\nsignals in HPV+ PSCC tumors compared to HPV- tumors, as shown in Revised Figure 4J-K.\nThis finding suggests that the interaction between IgA and PIGR is enhanced in HPV+ PSCC,\npotentially contributing to the differential immune response observed between HPV+ and HPV-\ntumors. We have included these new data in the revised manuscript to provide a more\ncomprehensive understanding of IgA-PIGR interactions in the tumor microenvironment.\nRevised Figure 4J-K. Colocalization of IgA and PIGR in PSCC tumor sections\n(J) Representative multiplex immunofluorescence images showing the distribution of PIGR and\nIgA in HPV- and HPV+ PSCC tumor. Colocalization of PIGR (green) and IgA (red), along with\nDAPI+ nuclei (blue), is indicated. Scale: 300 μm. (K) Quantification of IgA-PIGR colocalization\nevents per tumor area in HPV− (n=20) and HPV+ (n=20) PSCC tumors.\nFigure 6\nQ18: Figure 6N is inconsistent with the description in the text. The number of patients with double\nhigher IgA and PIGR expression and double lower IgA and PIGR expression is inconsistent in\nFigures 6N and 6O.\nResponse: We thank the reviewer for pointing out this inconsistency between Figure 6N and its\ndescription in the text. We have reviewed the figures and identified a writing error in the original\nmanuscript that caused this discrepancy. We have corrected this error in the revised manuscript,\nand the numbers now match appropriately.\nThe manuscript entitled: “A B cells-IgA-epithelial axis enhances antitumor response and improves\noutcome in HPV-associated penile squamous cell carcinoma” by Tao et al. aims to elucidate the\nmechanisms behind the improved clinical outcomes of HPV+ vs. HPVPSCC patients. The study\nanalyzes samples from a cohort of 90 PSCC patients using scRNA-Seq, scBCRSeq and bulk\nBCR-Seq complemented by IHC and validation in xenograft tumor models. The finding reveal that\nHPV+ PSCC exhibits elevated levels of IgA and polymeric Ig receptor corelating with improved\npatient survival compared to HPV- patients. Additionally, HPV+ PSCC patients displayed an\nenrichment of B cell and TLS in the TME. The report is well written, employs robust methodologies\nand presents convincing data. This study is important and opens new research avenues towards\nthe development of new therapeutic strategies in cancer treatment where the concept of\nharnessing both the T and B cell arms can enhance the anti-tumor impact. Nevertheless, there\nare some moderate issues that the authors should address as follows:\nResponse: We sincerely thank the reviewer for the valuable comments and for highlighting the\nsignificance of our work. We have carefully considered the issues raised and have addressed\nthem in the revised manuscript.\nQ1: B cell markers: There is an inconsistency in the definition of B cell markers. For example, in\nLine 158: the authors define B cells by the markers MS4A1+ and CD79A+ and plasma cells by J\nchain+ and IgG1/IgA1; Line 213: TIL-B are defined as MS4A1+ and plasma SDC1+; Figure 2C\nand Figure 2F: alternate between CD19 used as a marker in Figure 2C and CD20 used in Figure\n2F. How will this alternation affect the insights as these markers are expressed differentially during\nB cell development stages. Please clarify\nResponse: We thank the reviewer for highlighting the need to clarify our use of different B cell\nmarkers throughout the manuscript. Below, we explain the rationale for selecting specific markers\nfor B cells and plasma cells in each experimental approach and how these choices do not affect\nthe overall insights of the study:\n(1) Definition of B cells and plasma cells in scRNA-seq: In the revised scRNA-seq analysis,\nwe have standardized the definition, using MS4A1 (CD20) for B cells and SDC1 (CD138) for\nplasma cells (Revised Figure 3A-B). This approach aligns with common practices for scRNA-\nseq studies of B cell populations and ensures clear distinctions of TIL-B cells at the single-cell\nlevel.\n(2) Use of markers in flow cytometry (FACS): In original Figure 2C, we used CD19 to identify\ntumor-infiltrating B cells (TIL-Bs) through flow cytometry. CD19 is a well-established marker\nexpressed on most B cell subsets, including naïve, activated, and some plasma cells, but is\ntypically lost in fully differentiated plasma cells. We selected CD19 for FACS to broadly capture a\nwide range of B cells at various stages of development, ensuring a comprehensive evaluation of\nB cells, including subsets that may not express CD20.\n(3) Marker selection for IHC and TLS evaluation: In original Figure 2F, we used CD20 (MS4A1)\nfor immunohistochemistry (IHC) staining to assess tertiary lymphoid structures (TLSs). CD20 is a\nmature B cell marker commonly used for histological studies, especially in the context of tissue\nsections. This contrasts with CD19, which is less frequently used for IHC. The use of CD20\ninstead of CD19 in IHC was intentional, as CD20 more specifically identifies mature B cells within\nTLSs and is appropriate for evaluating their presence in tissue.\npopulations or structures in the relevant experimental context. These choices were\nmethodologically appropriate and do not affect the overall insights of the study. We have now\nclarified these marker definitions in the revised manuscript to ensure consistency and avoid\nconfusion.\nQ2: Figure 1E – it is not clear if the PN and PT are from the same patient - for example, PN5 is\naligned with PT4 - does this mean the PN4 was not available? also PN6 is aligned with PT5. This\nis confusing in a way that it hard to understand if the tissues are from the same patient – I\nrecommend to align the graph according to the patient they came from (leave an empty gap for\nPN4).\nResponse: We thank the reviewer for pointing out the confusion regarding the alignment of PN\nand PT samples in original Figure 1E.\nTo clarify, we collected fresh tumor specimens (PT1-PT6) and paired normal adjacent tissues\n(NAT) (PN1-PN6) from 6 treatment-naïve PSCC patients (P1-P6) for 3' scRNA-seq analysis in\nDiscovery Cohort 1. However, the 3' scRNA-seq library construction for the PN4 sample failed\ndue to an excessive number of dead cells, which is why PN4 data is missing.\nWe have updated the current alignment of samples in Revised Figure 1D to indicate that the\ncorresponding normal tissue sample from patient P4 was not available. This adjustment helps\nclarify that each PN and PT pair, where available, corresponds to the same patient.\nRevised Figure 1D-1E. The single-cell landscape of the tumor microenvironment in PSCC\n(D) The relative abundance of four cell components: immune cells, epithelial cells, endothelial cells and\ncancer-associated fibroblasts (CAF) in tumor tissues (n=6) or NATs (n=5) from individual PSCC patients,\nexcluding PN4. (E) Bar graph showing the relative proportions of the four cell components in tumor versus\nNATs.\nQ3: Although the overall distribution of cells indicates the dominance of immune cells in the TME\n(Figure 1D), looking at the distribution at the individual patient level there are some exceptions.\nFor example, PT1 and PN5. Please suggest a possible explanation for these differential\nsignatures.\nimmune cell signatures observed at the individual patient level in original Figure 1E. Upon\nreviewing the data, we noted the following observations:\nFor patient P1, both the tumor tissue (PT1) and the adjacent normal tissue (PN1) exhibited\nrelatively low immune cell infiltration. For patient P5, both the tumor tissue (PT5) and the adjacent\nnormal tissue (PN5) showed higher percentages of immune cells. We suggest that the differential\nimmune signatures in these patients may be explained by several contributing factors:\n(1) Patient-specific factors: Age and tumor size. P5 is an older patient with the largest tumor\nsize in our cohort. Age-related immune system changes, such as chronic inflammation or immune\ndysregulation, could contribute to increased immune cell infiltration in both the tumor and\nsurrounding normal tissues. Larger tumors also tend to induce more robust immune responses\ndue to greater antigen exposure and inflammation.\n(2) Tumor characteristics and immune evasion: (A) Low immune infiltration in P1. The reduced\nimmune cell infiltration in PT1 and PN1 may result from tumor immune evasion, such as PD-L1\noverexpression or secretion of immunosuppressive cytokines, suppressing immune response in\nboth tumor and adjacent tissues. (B) Higher tumor cell density in P1. P1 showed a higher\nproportion of malignant epithelial and endothelial cells, indicating a denser tumor\nmicroenvironment that could limit immune cell infiltration, possibly reflecting a more aggressive\nor invasive tumor state.\n(3) Heterogeneity in immune response: (A) High immune infiltration in P5. The high immune\ncell infiltration in both PT5 and PN5 could indicate a uniformly active immune response in this\npatient. The higher immune cell presence in both tumor and normal tissue could be a sign of an\nongoing immune surveillance mechanism that is active even in the adjacent normal tissue,\npotentially driven by chronic inflammation or a strong anti-tumor immune response. (B) Localized\nimmune suppression in P1. For P1, immune suppression may be more localized, leading to less\nimmune infiltration in both tumor and adjacent tissues, which may represent an immunologically\n\"cold\" tumor microenvironment.\n(4) Technical factors and sampling variability: The location of tissue sampling can also\ninfluence the observed immune cell composition. Small variations in the sampled regions of both\nthe tumor and the adjacent normal tissue could lead to different immune cell densities due to the\ninherent heterogeneity within and around the tumor microenvironment.\nRevised Figure S1A. Clinical information for PSCC patients that subjected to scRNA-seq\n(A) Clinical data for the PSCC patients is summarized in the table, detailing age, tumor size,\ntumor stage, in-situ hybridization (ISH) results for HPV16/18 DNA, and Immunohistochemistry\n(IHC) staining for p16. P1-P6 represent 6 treatment-naïve PSCC patients from whom we\ncollected both tumor specimens and paired normal adjacent tissue (NAT) for scRNA-Seq\n41\n(discovery cohort 1). P7-P10 represent 4 treatment-naïve PSCC patients from whom CD45+\nimmune cells were isolated from tumor specimens for scRNA-seq and paired B cell receptor\nsequencing (discovery cohort 2).\nQ4: Figure 1J – it is hard to differentiate between the clusters as the colors seem similar. Try to\nprovide color with high differentiative index. Also, maybe add a small arrow to highlight the\nnonimmune cells (cluster 22).\nResponse: We thank the reviewer for this suggestion. In our revised manuscript, we reanalyzed\nour scRNA-seq dataset using updated clustering parameters (20 principle components and a\nresolution of 0.5) to capture more subtle variation and avoid over-clustering, which led to the\nidentification of 20 clusters. We have adjusted the colors used for different clusters to ensure a\nhigh differentiative index. Additionally, as suggested, we have added clusters label to highlight\nthe non-immune cells. These changes have been made in the Revised Revised Figure 1I-1J to\nenhance the figure’s readability and improve the overall presentation of our data.\nRevised Figure 1. The single-cell landscape of the tumor microenvironment in PSCC\n(I) UMAP plot of 15,566 immune cells from the scRNA-seq datasets of discovery cohort 1,\ndisplaying 20 distinct immune clusters in PSCC tumor tissues and normal adjacent tissues (NATs).\n(J) Violin plot illustrating the normalized expression levels of selected marker genes for each\nidentified immune cluster.\nQ5: No pointer to figure 2C when mentioning the FACS.\nResponse: We thank the reviewer for pointing out this oversight. We have updated the FACS\ndata to improve the clarity in Revised Figure 2C.\nQ6: Figure 2D – Relative proportion out of CD45+? it so, it should be clearfied in the figure caption.\nPlease check all pointers to figure 2 panel as they seem off\nResponse: We thank the reviewer for this question.\nIn response, we confirm that the relative proportions displayed in the original Figure 2D were\nindeed calculated as a subset of CD45+ immune cells. To ensure this is clear, we have revised\nthe figure legend to explicitly state that the proportions are gated from CD45+ immune cells,\nclarifying that these represent immune populations within the broader CD45+ compartment\n(Revised Figure 2D).\nWe have thoroughly reviewed and corrected all figure pointers throughout the manuscript to\nensure clarity.\nRevised Figure 2C-2D. Tumor-infiltrating B cells and T cells in HPV+ and HPV- PSCC\n(C) Gating strategy used for FACS analysis of tumor-infiltrating B cells (TIL-Bs) and\ntumor-infiltrating T (TIL-Ts) cells from fresh tumor specimens of PSCC. CD45+ immune\ncells were selected, and TIL-Bs and TIL-Ts were further gated for detailed analysis. (D)\nHistograms showing the relative proportions of immune cells, including TIL-Bs, TIL-Ts,\nCD4+ TIL-Ts and CD8+ TIL-Ts as a fraction of the CD45+ population in fresh tumor\nspecimens from PSCC patients (HPV−, n=11; HPV+, n=9).\nQ7: Figure 2E – should described what type of statistical test was employed\nResponse: We thank the reviewer for this suggestion.\nIn response to the request for clarification, we confirm that the statistical analysis in original Figure\n2E, which presents survival curves comparing HPV+ and HPV- groups, was conducted using a\ntwo-sided log-rank (Mantel-Cox) test to assess statistical significance. We have now revised the\nfigure legend to clearly indicate the specific statistical test used.\nQ8: Figure 2G: how was the number of TLS counted? can the author comment on how will the\nsize of the tissue sample collected affect the number of TLS observed? If you count total TLS no\nin the context of size (as done in figure 2J) this probably will affect the number of TLSs identified\n(this is true for Figure 2G and 2I).\nnormalization of tertiary lymphoid structures (TLSs) in our analysis.\nTo accurately assess the number of TLSs, we followed a standardized procedure in which we\ncounted TLSs in immunohistochemistry (IHC) sections and normalized these counts by the tumor\narea. Normalizing by tumor area ensures that variations in tissue sample size do not introduce\nbias into the results. Only TLSs consisting of more than 50 cells were included in the analysis, as\nthis threshold ensures we are identifying significant TLS structures. This approach aligns with\nestablished methodologies for TLS quantification in previous studies9, 13, 14.\nIn the original manuscript, Figures 2G and 2I displayed the total number of TLSs without adjusting\nfor tumor area, while Figure 2J presented the number of TLSs normalized by tumor area. We\nagree with the reviewer’s concern that reporting the unnormalized TLS counts could introduce\nbias due to differences in tissue sample size. In response, we have revised the analysis to present\nonly the normalized number of TLSs in Revised Figures 2G and 2H, consistent with the approach\nused in Figure 2J. This change provides a more accurate representation of TLS density within the\ntumor microenvironment.\nRevised Figure 2G, 2H. Enrichment of Tertiary lymphoid structures in HPV+ PSCC tumor\n(G) Bar graph showing the number of tertiary lymphoid structures (TLSs) per square\ncentimeter in non-tumor (adjacent normal tissue) and intra-tumor regions of PSCC\npatients (n=37). (H) Bar graph comparing the number of TLSs per square centimeter in\nHPV− (n=52) and HPV+ (n=38) PSCC tumor sections. TLSs were counted in IHC-stained\nsections and normalized by tissue area to account for variations in sample size.\nQ9: Figure 2I: left graph, same question as for Figure 2G. It seems that when counting total there\nis not a significant difference between HPV- and HPV+ but when normalizing per area there is a\ndifference (Figure 2J). This should be also calrified in the text “…the number of TLSs and TIL-Bs\nin tumor tissues were significantly increased in HPV+ PSCC compared with HPV- controls”\nResponse: We thank the reviewer for this suggestion. In the latest version of the manuscript, we\nhave removed the figures presenting the total numbers and now show only the normalized\nnumbers of TLSs by tumor area. We have updated the text accordingly to reflect this more\nrigorous approach.\nResponse: We thank the reviewer for pointing out this oversight. We have updated the text in the\nmanuscript as below.\n‘We evaluated the spatial localization of B and T cells in tumor sections and observed the\npresence of tertiary lymphoid structures (TLSs) (Revised Figure 2F). TLSs are immune cell\naggregates consisting of B and T cells that form outside secondary lymphoid organs (SLOs) in\nresponse to chronic inflammation or infection. Previous studies have demonstrated that TLSs are\nassociated with improved patient survival in various human tumors, including melanoma, renal\ncell carcinoma, and sarcoma.’\nQ11: Line 232-237: Moreover, a comparison between the B cell composition in the TLS\nthemselves would be informative as well.\nResponse: We thank the reviewer for this insightful suggestion regarding the comparison of B\ncell composition within TLSs. In response to this, we conducted additional analyses to evaluate\nthe specific B cell subsets within TLSs using multiplex immunofluorescence (mIF) staining. This\nanalysis was performed on a random and unbiased cohort of 20 PSCC tumor sections, consisting\nof 10 HPV+ and 10 HPV- patients.\nTo gain a more comprehensive understanding of the B cell composition within TLSs, we selected\nkey markers for mIF staining: (1) CD3 to assess T cell populations and their interactions with B\ncells within TLSs, (2) CD20 to identify mature B cells, (3) CD38 to mark activated B cells and\nplasma cells, indicating their differentiation and activation status, and (4) BCL6 to identify germinal\ncenter B cells, which are crucial for antibody maturation and the overall functionality of TLSs.\nThese markers allowed us to assess the diversity and activation state of immune cells within TLSs\nin both HPV+ and HPV- PSCC tumors.\nOur mIF staining produced clear and reliable signals for CD3, CD20, CD38. We observed well-\nstructured TLS, with CD3 and CD20 staining clearly outlining their organization, more distinctly\nthan our previous IHC staining. Notably, PNAd+ cells were detected within TLS, though the number\nof PNAd+ cells was quite low, which is consistent with previously published findings (Immunity,\n2022). Unfortunately, the Bcl6 signal was relatively weak in our samples, making it difficult to\nclearly identify germinal centers. Therefore, we proceeded with analyzing only CD3, CD20, and\nCD38 staining.\nOur analysis showed that TLSs in HPV+ PSCC exhibited significantly higher numbers of CD20+,\nand CD38+ cells compared to TLSs in HPV- PSCC (normalized by tumor area to account for\nsample variability). These findings indicate that TLSs in HPV+ PSCC are more enriched in\nactivated B cells and plasma cells, reflecting a more robust immune response. The representative\nimages and quantification of the findings are presented in Revised Figure 2I-2J.\n(I) Representative multiplex immunofluorescence images showing the distribution of four key cell\ntypes within TLSs in PSCC tumor section. Markers include CD3+ T cells, CD20+ B cells, CD38+\nplasma cells, BCL6+ germinal center B cells, and DAPI+ nuclei. Scale bar: 300 μm. (J)\nQuantification of the numbers of CD3+, CD20+, and CD38+ cells within TLS in HPV− (n=10) and\nHPV+ (n=10) PSCC tumor samples. Data are presented as the number of cells within TLS,\nnormalized by tumor area for comparative analysis.\nQ12: Line 961 (and figure 3C-D): The way a B cell clone is defined is not clear enough. The way\nclones are defined can have a substantial impact on the biological insights. It seems that they\nemployed a simplified strategy to clone the B cells by counting the number of unique VH and\ndefining clonality by those that exhibit more the one BCR. This should be revised according to a\nmore acceptable clonotyping strategy. For example, same V gene+ same J gene + identical\nCDRH3. Noteworthy, there are many clonotyping strategies reported in literature, the authors\nshould use one of these strategies or proposed a new strategy which has a biological rational.\nResponse: We sincerely appreciate the reviewer’s suggestion and have taken steps to ensure\nthat the manuscript now provides a detailed, robust explanation of our methodology for defining\nB cell clones in our scBCR-seq analysis.\nIn our study, we analyzed B cell clonality using the scRepertoire package (version 1.3.5), which\nis widely used for this purpose7. Clonotypes were assigned only for cells where high-quality paired\nheavy and light chains were sequenced. This is essential for accurately defining B cell clones and\navoiding misclassification. To ensure a biologically meaningful and stringent clonotyping strategy,\nwe utilized the CTstrict function within scRepertoire, which defines clones based on the following\ncriteria: (1) V(D)J gene usage. Clones were defined by identifying B cells that share the same V\ngene and J gene usage in their immunoglobulin sequences. This method ensures that clonally\nrelated B cells have undergone the same recombination events during B cell development. (2)\nCDR3 sequence identity. The CDR3 (complementarity-determining region 3), the most variable\nand functionally significant region of the B cell receptor, was compared across cells. To define\nclonally related cells, we used a >85% normalized hamming distance for the nucleotide sequence\nof the CDR3 region. This threshold ensures that minor variations due to sequencing errors or\nsomatic hypermutation are accounted for, while maintaining high stringency for identifying clonal\nrelationships.\nBy employing this combination of V(D)J gene usage and CDR3 sequence identity, we ensured\nthat our clonotyping strategy was rigorous and aligned with best practices in the field. This method\ngoes beyond simply counting unique VH sequences, as it also includes key determinants such\nas J gene usage and CDR3 sequence similarity, which are critical for defining B cell clones with\nhigh confidence. Our methodology is consistent with recent studies that have employed similar\nstrategies to assess B cell clonality using scBCR-seq8. We have updated the Methods section\nin our manuscript to more clearly articulate this approach.\na reference that describes the usage of EGR1 as a marker for plasma cells?\nResponse: We thank the reviewer for this question regarding our cell classification strategy. We\nhave reanalyzed tumor-infiltrating B cells (TIL-Bs) in scRNA-seq data.\nAligning with the reviewers’ suggestions, we included CD138 (SDC1) as a well-established\nmarker for plasma cells, which typically lose expression of CD19 and CD20 as they mature. So,\nplasma cells are defined by the expression of CD138 (SDC1) along with high levels of CD38.\nBesides, to differentiate plasmablasts, we used a combination of CD20 (MS4A1), CD38, and\nMKI67 markers. Plasmablasts are defined as CD38+ but lack CD138 expression and are\nCD20loMKI67+. Moreover, we used additional markers like CD24, CD27, and EGR1 to distinguish\nbetween different memory B cell subsets, rather than relying solely on immunoglobulin heavy\nchain genes (PMID: 35393541, Panel C)24.\nRegarding EGR1, we did not use this marker to define plasma cells. Instead, we identified two\ndistinct clusters of memory B cells based on the expression of EGR1: Memory B cells-1 (EGR1-)\nand Memory B cells-2 (EGR1+). EGR1 has been reported to play a non-redundant role in the\ndifferentiation of B cells into plasma cells (PMID: 26140048)25, indicating its involvement in the\nearly stages of B cell differentiation rather than in terminally differentiated plasma cells. This\ndifferentiation into memory B cells based on EGR1 expression is also supported by studies that\nshow EGR1 expression in antigen-specific activated B cell subtypes committed to becoming\nmemory B cells (scRNA-seq of oesophageal squamous cell carcinoma, PMID: 36650114) 26.\nPanel D\nBased on this reanalysis, we have categorized TIL-Bs into six clusters:\nnaïve B cells (C1, MS4A1+IGHD+IGHM+)\nEGR1- or EGR1+ memory B cells (C2-C3, MS4A1+IGHD-IGHM-CD24+CD27+)\ngerminal center B cells (C4, MS4A1+MME+)\nplasmablast (C5, MS4A1loCD38+SDC1-MKI67+)\nplasma cells (C6, MS4A1-CD38+SDC1+)\nWe have made these adjustments to our analysis and updated in Revised Figure 3A-3B\naccordingly to ensure that our classification of B cell subsets and plasma cells is accurate and\nreflects current best practices in the field.\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing), illustrating six clusters. The clusters are coloured by cell type (left), individual patient\n(middle) and HPV status (right). (B) Violin plots showing the normalized expression of selected\nmarker genes corresponding to each TIL-Bs cluster. (C) Stages of B cell activation, clonal\nexpansion, somatic hypermutation, affinity selection, class switch recombination and\ndifferentiation to plasma cells (PCs). (D) Heatmap showed representative genes of B cells in each\nsubtype in oesophageal squamous cell carcinoma (ESCC)\nQ14: Line 260: B cell clonality - it not clear to me how B cell clonality is determined? is B cell\nclonality determined if there are more than one unique heavy chain variable region sequencing in\na cluster? please clarify. Moreover, how is the B clonal diversity determined?\nResponse: We thank the reviewer for this insightful question, and we apologize for any confusion\ncaused. We would like to provide a clearer explanation regarding how B cell clonality and clonal\ndiversity were determined in our analysis.\n1. Defining B cell clonality: As detailed in our response to Q12, we assessed B cell clonality\nusing the CTstrict function from the scRepertoire package (version 1.3.5), a widely used tool for\nclonotype analysis. B cell clones were identified based on two key criteria: (1) V(D)J gene usage:\nWe grouped B cells that shared the same V gene and J gene usage in their immunoglobulin\nsequences. This ensures that the B cells being classified into clones originated from the same\nrecombination events. (2) CDR3 sequence identity: We used a >85% normalized hamming\ndistance to define clonal relationships, which accounts for minor variations due to somatic\nhypermutation or sequencing errors while maintaining stringency. This combination of V(D)J gene\nusage and CDR3 sequence similarity ensures a rigorous identification of B cell clones and aligns\nwith best practices in the field8.\n2.Clonal diversity: Clonal diversity refers to the variety and distribution of distinct B cell clones\nwithin a given population, providing insight into the complexity of the B cell repertoire. To assess\ndiversity, we used the Shannon diversity index, a widely accepted metric for measuring diversity\nin biological systems. The Shannon diversity index accounts for both the number of unique BCR\nsequences (i.e., the richness of the repertoire) and the evenness of their distribution across\ndifferent clones. We calculated this index using the clonalDiversity function from the scRepertoire\npackage, which allowed us to evaluate the diversity of BCRs within each cluster and gain insights\ninto the immune complexity and the extent of clonal expansion.\nmanuscript to further clarify the details of our clonotyping strategy and the calculation of clonal\ndiversity.\nQ15: Line 304-306: The author wrote: “Taken together, these data suggest that HPV+ PSCC local\ntumor tissues exhibited higher levels of IgA, which correlated with improved survival in PSCC\npatients.”. From the way this sentence is phrased it seems that the authors are no convinced that\ntumor tissues exhibited higher levels of IgA?\nResponse: We thank the reviewer for this observation. We apologize if the original phrasing\ncaused any confusion. We are confident in our finding that HPV+ PSCC tumor tissues exhibit\nhigher levels of IgA compared to HPV- tumors.\nTo clarify this point, we have revised the statement in the manuscript to accurately reflect our\nfindings: These data demonstrate that HPV+ PSCC local tumor tissues exhibit higher levels of IgA\ncompared to HPV- tumors, and this increase is correlated with improved survival in PSCC\npatients.\nQ16: Line 904: There is no mention of steps taken to remove ambient RNA, which may be\nimportant to ensure the accuracy and quality of the data. Ambient RNA can result from cell lysis\nduring sample preparation and can significantly contaminate the transcriptomic profiles of\ncaptured cells (create background noise). Including a method for depleting this RNA can improve\nthe reliability of the sequencing results. Please clarify and explain why ambient RNA was not\nremoved.\nResponse: We thank the reviewer for this question regarding the removal of ambient RNA to\nensure the accuracy and quality of our scRNA-seq data.\nWe acknowledge that ambient RNA, which can result from cell lysis during sample preparation,\nhas the potential to contaminate transcriptomic profiles and introduce background noise. To\naddress this concern, we implemented several steps to evaluate and mitigate the impact of\nambient RNA, ensuring the reliability of our sequencing results:\n(1) Initial quality assessment. Before proceeding with downstream analysis, we conducted\nmultiple quality control assessments: (A) CellRanger output alerts: we closely examined the\nCellRanger output for any \"Low Fraction Reads in Cells\" alerts, which indicate potential ambient\nRNA contamination. No such alerts were found in our data, suggesting minimal contamination.\n(B) Barcode rank plots: we reviewed barcode rank plots for a characteristic ‘steep cliff,’ which\ntypically signals a clear distinction between high-quality cells and background noise. Our data\nexhibited this expected pattern, further indicating the absence of significant contamination. (C)\nMitochondrial gene enrichment: we analyzed the marker genes for each cluster to detect any\nunusual enrichment of mitochondrial genes, which could be indicative of ambient RNA\ncontamination. Our analysis showed no significant mitochondrial gene enrichment.\n(2) Reanalysis using SoupX. Although the initial assessments indicated minimal risk of ambient\nRNA contamination, we took an additional step by reanalyzing the data using SoupX27 (PMID:\n33367645), a tool specifically designed to identify and remove ambient RNA contamination.\nSoupX confirmed that the levels of ambient RNA were negligible and did not significantly impact\nthe transcriptomic profiles of our captured cells. No additional cells were flagged for removal\nbeyond those filtered during our initial quality control steps (Figure below).\nthat ambient RNA contamination was not a significant factor in our dataset. The steps we took to\nevaluate ambient RNA ensured that the transcriptomic data used in our analyses is of high quality\nand accuracy. We have updated the manuscript to include this explanation and to clarify our\nmethodology in addressing potential ambient RNA contamination.\nQ17: General comment: The figure captions require enhancement to better guide the reader\nthrough the data presented in the figures. Currently, several captions lack sufficient detail, which\nhinders comprehension. For instance, Figure 2D is captioned: \"Histograms showing the relative\nproportions of TIL-Bs and T cells in 20 fresh tumor specimens of PSCC (HPV−, n=11; HPV+,\nn=9).\" However, it remains unclear what the baseline cell population is, among other details.\nPlease ensure that each caption contains all necessary information to allow for independent\ninterpretation of the figures without referring back to the main text.\nResponse: We thank the reviewer for this valuable feedback regarding the figure captions.\nWe recognize the importance of providing detailed and comprehensive figure captions to guide\nreaders through the data presented in the figures. To address this, we have carefully revised all\nfigure captions to ensure they include all necessary information for independent interpretation,\nwithout needing to refer back to the main text.\nFor example, the updated caption for Figure 2D now reads: Histograms showing the relative\nproportions of immune cells, including TIL-Bs, TIL-Ts, CD4+ TIL-Ts and CD8+ TIL-Ts as a\nfraction of the CD45+ population in 20 fresh tumor specimens from PSCC patients (HPV−,\nn=11; HPV+, n=9).\nWe have applied similar improvements to other figure captions to ensure they provide sufficient\ndetails for understanding the figures independently. We appreciate the reviewer’s suggestion and\nbelieve the revise will improve the overall readability and comprehension of our manuscript.\nUsing a cohort of a cohort of 90 patients with penile cancer, Tao et al have characterized the\nimmune-environment of this disease, comparing HPV infection-dependent differences. The\nauthors found that HPV+ penile cancers are heavily enriched in intra-tumoral B cells, plasma cells,\nand TLSs, with abundant IgA secretion (predictive of outcome) and expression of polymeric IgA\nreceptor (PIGR) in tumor cells. Higher levels of IgA and PIGR expression are associated with\nbetter outcome, whereas the effects of IgA:PIGR interactions are characterized. The manuscript\nrecapitulates many of the concepts recently reported for human ovarian and endometrial cancer\n(and in NSCLC to a lesser extent), although those manuscripts are blatantly ignored (e.g., PMID:\n33536615, 34949671, 37909039). These data are nevertheless relevant because: 1) Little is\nknown about the ME of penile tumors; much less in a cohort of this size; 2) the results represent\nthe first independent confirmation of mechanisms reported in other tumors by a single group\n(although ignored in the references); 3) there are significant differences between HPV+ and HPV-\ntumors for B cell infiltration and IgA production, but not for other immune determinants; 4) the\nauthors made extensive use of single cell OMICS, with interesting differences across HPV-driven\ntumors; and 5) the unusual density of TLS in HPV+ tumors is interesting and could illustrate the\nimmunogenicity of these tumors. The strength of this manuscript is in the human data. The mouse\ndata are questionable and need to be confirmed, or ignored. In addition to adding the references\npreviously ignored, the authors must address the following major comments:\nResponse: We sincerely thank the reviewer for the thoughtful and comprehensive feedback.\nWe are grateful for the recognition of the value of our human data, especially given the rarity of\npenile cancer and the size of our cohort. Besides, we acknowledge the oversight in not citing the\nrecent relevant studies (e.g., PMID: 33536615, 34949671, 37909039) that have explored similar\nimmune mechanisms in other tumor types. We have now updated the manuscript to incorporate\nthese important references, providing a more thorough discussion of the parallels and differences\nbetween our findings in PSCC and the existing data in other cancers. Moreover, we have carefully\nre-examined the results from our mouse models in response to the reviewer’s concerns. Where\nneeded, we have reinforced the experimental validation to improve the robustness of the data or\nlimited our conclusions to avoid overstating their implications, ensuring that the focus remains on\nthe solid human data presented.\nQ1: The distinction of B cells and plasma cells in Figure 1 seems to be based on expression of\nIGH1 and IGHA1 to define the latter. Independently of the possibility that plasma cells could\nproduce other isotypes (e.g., IGHA2 or IGG4), isotype switched B cells can also express those\ngenes. CD138 expression (as shown at the protein level in Figure 3) would be a definitory marker\nof plasma cells, while CD138- plasmoblasts could be defined as CD19+CD20-CD38+CD27-high\ncells. in any case, better markers are needed to establish these categories, because >10% of\nplasma cells in tumor beds (Fig1L) is unusually high.\nResponse:\nWe thank the reviewer for this insightful question. We understand the importance of using precise\nmarkers to differentiate between various B cell subsets, particularly in distinguishing plasma cells\nfrom other B cells. In response to the reviewer’s suggestion, we have reanalyzed the tumor-\ninfiltrating B cells (TIL-Bs) using a more refined set of markers.\nmarker for plasma cells, which typically lose expression of CD19 and CD20 as they mature. So,\nplasma cells are defined by the expression of CD138 (SDC1) along with high levels of CD38. This\nensures that we are accurately identifying mature plasma cells within the tumor microenvironment.\nBesides, to differentiate plasmablasts, we used a combination of CD20 (MS4A1), CD38, and\nMKI67 markers. Plasmablasts are defined as CD38+ but lack CD138 expression and are\nCD20loMKI67+. Moreover, we used additional markers like CD24, CD27, and EGR1 to distinguish\nbetween different memory B cell subsets, rather than relying solely on immunoglobulin heavy\nchain genes.\nBased on this reanalysis, we have categorized TIL-Bs into six clusters:\nnaïve B cells (C1, MS4A1+IGHD+IGHM+)\nEGR1- or EGR1+ memory B cells (C2-C3, MS4A1+IGHD-IGHM-CD24+CD27+)\ngerminal center B cells (C4, MS4A1+MME+)\nplasmablast (C5, MS4A1loCD38+SDC1-MKI67+)\nplasma cells (C6, MS4A1-CD38+SDC1+)\nWe have made these adjustments to our analysis and updated in Revised Figure 3A-3B\naccordingly to ensure that our classification of B cell subsets and plasma cells is accurate and\nreflects current best practices in the field. The updated data provide a more precise\ncharacterization of the B cell composition within the tumor microenvironment, addressing the\nconcern of unusually high plasma cell percentages.\nRevised Figure 3A-3B. scRNA-seq analysis of TIL-Bs in PSCC\n(A) UMAP plot of 2,362 TIL-Bs from discovery cohort 2 (5’ scRNA-seq paired with scBCR-\nsequencing), illustrating six clusters. The clusters are coloured by cell type (left), individual patient\n(middle) and HPV status (right). (B) Violin plots showing the normalized expression of selected\nmarker genes corresponding to each TIL-Bs cluster.\nQ2: The structures shown in Figure 2 are indeed very strongly suggestive of TLS, but the\npresence of other components such as PNAd+ high endothelial venules and, ideally, CD21\nfollicular dendritic cells, would confirm this definition.\nResponse:\nendothelial venules and CD21+ follicular dendritic cells (FDCs), to confirm the presence of mature\nTLS. To address this, we conducted additional multiplex immunofluorescence (mIF) staining on\nPNAd, CD21, DAPI and (2) CD3, CD20, CD38, Bcl6, DAPI, to evaluate tertiary lymphoid\nstructures (TLS) and the composition of B cells within these structures.\nIn this analysis, PNAd+ high endothelial venules and CD21+ FDCs were specifically included, as\nrecommended, to verify the presence of mature TLS. Additionally, CD38 was used to identify\nactivated B cells and plasma cells, and BCL6 to mark germinal center B cells, providing further\ninsight into the B cell composition within TLS.\nOur mIF staining produced clear and reliable signals for CD3, CD20, CD38, and PNAd. We\nobserved well-structured TLS, with CD3 and CD20 staining clearly outlining their organization,\nmore distinctly than our previous IHC staining. Notably, PNAd+ cells were detected within TLS,\nthough the number of PNAd+ cells was quite low, which is consistent with previously published\nfindings (Immunity, 2022).\nUnfortunately, we faced technical challenges with the CD21 staining. Despite testing two different\nCD21 antibodies—mouse-anti-human CD21 (Servicebio, GB14031-50) and mouse-anti-human\nCD21 (ZsBio, ZM0040)—we were unable to obtain reliable staining signals for CD21+ FDCs.\nAdditionally, the Bcl6 signal was relatively weak in our samples, making it difficult to clearly identify\ngerminal centers. Therefore, we proceeded with analyzing only CD3, CD20, and CD38 staining.\nWe have updated the manuscript in Revised Figure 2I, 2J to reflect these findings, we\nacknowledge the limitations in our ability to reliably stain CD21+ FDCs and are committed to\nfurther optimizing the staining protocols for future studies.\nRepresentative multiplex immunofluorescence images showing the distribution of key cell types\nwithin TLSs in PSCC tumor section. Two marker sets include: (Top Right) CD3+ T cells, CD20+\nB cells, PNAd+ cells, and DAPI+ nuclei, scale:200μm; (I) CD3+ T cells, CD20+ B cells, CD38+\nplasma cells, BCL6+ germinal center B cells, and DAPI+ nuclei. Scale bar: 300 μm. (J)\nQuantification of the numbers of CD3+, CD20+, and CD38+ cells within TLS in HPV− (n=10) and\nHPV+ (n=10) PSCC tumor samples. Data are presented as the mean number of cells within TLS,\nnormalized by tumor area for comparative analysis.\nQ3: The data shown in Fig5D-G are puzzling. It has been reported that PIGR provides an\noncogenic signal in at least hepatocarcinoma by activating ribosomes through a poorly\ncharacterized mechanism (PMID: 33390805). If expression of PIGR is always deleterious for the\ntumor, even in the absence of IgA, why tumor cells do not evolve to lose PIGR?\nResponse: We thank the reviewer for raising this insightful question regarding the role of PIGR\nin cancer. We recognize that PIGR has been reported to play diverse roles in different cancer\ntypes, exhibiting both pro-tumorigenic and anti-tumorigenic functions depending on the context.\nFor example, PIGR has been associated with pro-tumor functions in hepatocarcinoma28, 29. In\ncontrast, PIGR has been shown to exhibit anti-tumor effects in ovarian cancer30 and lung cancer31.\nSpecifically, in A549 lung cancer cells and 16-HBE cells, overexpression of PIGR inhibited cell\nproliferation in vitro31, which aligns with our findings in original Figure 5D-5G. These data suggest\nthat PIGR’s role in tumor progression is highly context-dependent and may vary significantly\nacross cancer types.\nIn our study, HPV- PSCC tumors indeed exhibited low PIGR expression, which align with the\nreviewer’s point that tumor evolution can favor the downregulation of PIGR to promote tumor\nprogression. However, in HPV+ PSCC tumors, PIGR expression was significantly elevated, driven\nby pro-inflammatory cytokines such as IFN-γ, which was also upregulated in HPV+ tumors\ncompared to HPV- tumors (Revised Figures 4M-4P). This suggests that PIGR's expression and\nfunction in HPV+ PSCC are shaped by the inflammatory tumor microenvironment, where chronic\ninflammation plays a prominent role.\nThe retention of PIGR expression in HPV+ PSCCs, despite its anti-tumor roles, may serve specific\npurposes, such as: (1) Modulating the immune microenvironment: PIGR may help regulate the\nlocal immune microenvironment, ensuring a controlled immune response rather than allowing an\nuncontrolled, aggressive immune reaction that kill tumor.(2) Maintaining epithelial integrity: PIGR\nmay help preserve epithelial integrity and barrier function, which can create a microenvironment\n54\nthat allows the tumor to evade more destructive immune responses in the context of chronic\ninflammation.\nThus, the role of PIGR in cancer is complex and likely influenced by both the tumor’s intrinsic\ncharacteristics and the external pressures of its microenvironment. Whether PIGR expression is\nretained or lost depends on the tumor’s need to balance immune regulation, inflammation, and\ntissue integrity. We have updated the manuscript to incorporate this nuanced view, emphasizing\nthe context-dependent nature of PIGR’s role across different cancer types, particularly in HPV+\nPSCC.\nQ4: For Fig5J-R, previous studies reported RNA-seq changes elicited by IgA transcytosis through\nendometrial (PMID: 34949671) and ovarian (PMID: 33536615) cancer cells. How do these results\ncompare to those analyses?\nResponse: We thank the reviewer for this insightful question and the opportunity to compare our\nfindings with previous studies on IgA transcytosis in cancer cells.\nTo address the reviewer’s query, we reanalyzed RNA-seq data from two relevant studies—on\nendometrial cancer (PMID: 34949671) and ovarian cancer (PMID: 33536615)—that investigated\nthe transcriptomic changes induced by IgA transcytosis.\nOur comparison revealed that primary PSCC cells treated with IgA showed significant enrichment\nin pathways related to the 'interferon-alpha response', 'interferon-gamma response' and\n'apoptosis' (Revised Figure 5O-Q), which is consistent with findings from both endometrial and\novarian cancer cells. This suggests a conserved mechanism across different cancer types, where\nIgA, mediated by PIGR transcytosis, triggers immune responses and promotes apoptosis,\npotentially contributing to anti-tumor effects.\nHowever, we also identified pathways that were uniquely downregulated in PSCC cells, such as\n‘cell cycle,’ ‘G2M checkpoint,’ and ‘E2F targets,’ which were not significantly affected in the\nendometrial and ovarian cancer studies. These PSCC-specific changes suggest that the IgA-\nPIGR axis may exert additional inhibitory effects on cell cycle regulation in PSCC, contributing to\nits anti-tumor activity in this specific cancer type.\nOverall, the transcriptional changes induced by IgA transcytosis highlight the role of IgA-PIGR\naxis in modulating the tumor immune microenvironment and enhancing tumor suppression\nthrough immune signaling and apoptosis. Our findings, in conjunction with existing research,\nsuggest that this pathway could serve as a therapeutic target in cancers with high PIGR\nexpression. We have updated the manuscript to reflect these comparative insights and to\nemphasize the broader relevance of our results.\n(Revised Figure 5O-Q, A, B) Gene set enrichment analysis (GSEA) illustrating the key pathways\nupregulated or downregulated in IgA-treated PSCC primary cells, (A) endometrial cancer cells,\nand (B) ovarian cancer cells. In PSCC cells, significant enrichment was observed in pathways\nrelated to the 'interferon-alpha response,' 'interferon-gamma response,' and 'apoptosis' within the\nIgA-PIGR group. Conversely, pathways related to the ‘cell cycle,’ ‘G2M checkpoint,’ and ‘E2F\ntargets’ were enriched in the IgA-control group. In both endometrial and ovarian cancer cells, IgA\ntreatment similarly induced the upregulation of immune response and apoptosis pathways,\nwithout notable alterations in cell cycle-related pathways.\noverexpressing Pigr in MC38, as opposed to CRISPR naturally expressed Pigr?; and 2) Do MC38\ntumors really induce the production of IgA responses? This needs to be shown and it would be\nweird, because rapidly progressing mouse models induce limited isotype switched responses.\nResponse: We thank the reviewer for these thoughtful questions.\n(1) Rationale for overexpressing PIGR in MC38 cells: In our study, we found that MC38 cells\nhave minimal endogenous PIGR expression. Since PIGR is significantly upregulated in HPV+\nPSCC and correlates with better overall survival, we chose to overexpress PIGR in MC38 cells to\nmodel its effects in a more differentiated tumor setting. A CRISPR knockout approach would not\nbe suitable here, as the baseline PIGR expression in MC38 cells is already very low, and knocking\nout PIGR would not provide meaningful insights into its role in tumor suppression.\n(2) IgA production in the MC38 tumor model: We acknowledge that subcutaneous MC38\ntumors typically have limited B cell infiltration, and, as the reviewer mentioned, rapidly progressing\ntumor models tend to have restricted isotype-switched immune responses. This is due to the short\ntime frame in which the immune system can activate B cells and undergo class-switch\nrecombination. However, we did observe some IgA within the MC38 tumors. Following\nadministration of anti-CD20 antibodies, we effectively depleted CD20+ B cells from the blood and\nspleen, which resulted in a reduced presence of IgA and IgG in tumor tissues (Revised Fig. 6E-\n6F). This suggests that IgA and IgG were present in the tumor microenvironment, although they\nare likely sourced from peripheral sites rather than produced locally by tumor-infiltrating B cells.\nThe IgA detected in the tumors may have originated from plasma cells in nearby lymphoid tissues,\nsuch as draining lymph nodes or mucosal-associated lymphoid tissues (MALT), which have a\nhigh capacity for IgA production. IgA could enter the tumor microenvironment from the systemic\ncirculation, where it is transcytosed into the tumor by PIGR-expressing epithelial or tumor cells.\nThis mechanism could explain the presence of IgA in tumors with limited B cell infiltration20, 21.\nWe recognize that the MC38 model has limitations in studying local IgA production. Indeed, our\ndata showed that the absence of IgA alone did not significantly affect MC38 tumor growth,\nindicating that IgA may not be the primary driver of tumor suppression in this model (Revised\nFigure 6H). However, our goal was to investigate the downstream effects of IgA, particularly its\ninteraction with PIGR. We observed enhanced anti-tumor effects with PIGR overexpression, and\nthese effects were significantly reduced in the absence of IgA. This suggests that the IgA-PIGR\naxis plays a critical role in modulating the tumor microenvironment and enhancing tumor\nsuppression (Revised Figure 6H).\nWe have updated the manuscript to clarify these points and provide a more detailed explanation\nof the experimental design and findings.\ngrowth\n(E, F) C57B6/J wild-type mice inoculated subcutaneously with 5´105 control- and PIGR-\noverexpressed MC38 cells. 250mg anti-CD20 antibodies (aCD20 mAb) per mouse were\nadministered at the time of MC38 tumor cell implantation and again at 10 days intervals\nto deplete B cells. n=4. ELISA quantification of the concentration of intra-tumor IgG (E)\nand intra-tumor IgA (F) in the four groups. (H) tumor growth curves in WT (IgA+/+) and IgA\ndeficient (IgA-/-) mice inoculated subcutaneously with 5´105 control- and PIGR-\noverexpressed MC38 cells. n=5. All data are shown as mean ± s.e.m. ns: not significant,\n*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001. Data are representative of two\nindependent experiments for (A-D) and (G-J).\nMINOR COMMENTS\nQ6: There is a typo in the title of Figure 5M\nResponse: We thank the reviewer for pointing out this oversight. We have corrected the typo in\nFigure 5M title.\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is\npart of the Nature Communications initiative to facilitate training in peer review and to provide\nappropriate recognition for Early Career Researchers who co-review manuscripts.\nReferences\n1. Emmanuel, A., Nettleton, J., Watkin, N. & Berney, D.M. The molecular pathogenesis of\npenile carcinoma-current developments and understanding. Virchows Arch 475, 397-405\n(2019).\n2. Ruffin, A.T. et al. B cell signatures and tertiary lymphoid structures contribute to outcome\nin head and neck squamous cell carcinoma. Nat Commun 12, 3349 (2021).\n3. Cillo, A.R. et al. Immune Landscape of Viral- and Carcinogen-Driven Head and Neck\nCancer. Immunity 52, 183-199 e189 (2020).\n4. Cubilla, A.L. et al. The World Health Organisation 2016 classification of penile\ncarcinomas: a review and update from the International Society of Urological Pathology\nexpert-driven recommendations. Histopathology 72, 893-904 (2018).\n5. Thomas, J. & Primeaux, T. Is p16 immunohistochemistry a more cost-effective method for\nidentification of human papilloma virus-associated head and neck squamous cell\ncarcinoma? Ann Diagn Pathol 16, 91-99 (2012).\n6. Steinestel, J. et al. The role of histologic subtype, p16(INK4a) expression, and presence of\nhuman papillomavirus DNA in penile squamous cell carcinoma. BMC Cancer 15, 220\n(2015).\n7. Borcherding, N., Bormann, N.L. & Kraus, G. scRepertoire: An R-based toolkit for single-\ncell immune receptor analysis. F1000Res 9, 47 (2020).\n8. Mathew, N.R. et al. Single-cell BCR and transcriptome analysis after influenza infection\nreveals spatiotemporal dynamics of antigen-specific B cells. Cell Rep 35, 109286 (2021).\n9. Vanhersecke, L. et al. Mature tertiary lymphoid structures predict immune checkpoint\ninhibitor efficacy in solid tumors independently of PD-L1 expression. Nat Cancer 2, 794-\n802 (2021).\n10. Cabrita, R. et al. Tertiary lymphoid structures improve immunotherapy and survival in\nmelanoma. Nature 577, 561-565 (2020).\n11. Petitprez, F. et al. B cells are associated with survival and immunotherapy response in\nsarcoma. Nature 577, 556-560 (2020).\n12. Helmink, B.A. et al. B cells and tertiary lymphoid structures promote immunotherapy\nresponse. Nature 577, 549-555 (2020).\n59\n13. Kinker, G.S. et al. B Cell Orchestration of Anti-tumor Immune Responses: A Matter of\nCell Localization and Communication. Front Cell Dev Biol 9, 678127 (2021).\n14. Meylan, M. et al. Tertiary lymphoid structures generate and propagate anti-tumor\nantibody-producing plasma cells in renal cell cancer. Immunity 55, 527-541 e525 (2022).\n15. Olesen, T.B. et al. Prevalence of human papillomavirus DNA and p16(INK4a) in penile\ncancer and penile intraepithelial neoplasia: a systematic review and meta-analysis. Lancet\nOncol 20, 145-158 (2019).\n16. Lohr, M. et al. The prognostic relevance of tumour-infiltrating plasma cells and\nimmunoglobulin kappa C indicates an important role of the humoral immune response in\nnon-small cell lung cancer. Cancer Lett 333, 222-228 (2013).\n17. Kroeger, D.R., Milne, K. & Nelson, B.H. Tumor-Infiltrating Plasma Cells Are Associated\nwith Tertiary Lymphoid Structures, Cytolytic T-Cell Responses, and Superior Prognosis\nin Ovarian Cancer. Clin Cancer Res 22, 3005-3015 (2016).\n18. Tan, X. et al. Molecular stratification by BCL2A1 and AIM2 provides additional\nprognostic value in penile squamous cell carcinoma. Theranostics 11, 1364-1376 (2021).\n19. Affara, N.I. et al. B cells regulate macrophage phenotype and response to chemotherapy in\nsquamous carcinomas. Cancer Cell 25, 809-821 (2014).\n20. Mora, J.R. & von Andrian, U.H. Differentiation and homing of IgA-secreting cells.\nMucosal Immunol 1, 96-109 (2008).\n21. Bemark, M. & Angeletti, D. Know your enemy or find your friend?-Induction of IgA at\nmucosal surfaces. Immunol Rev 303, 83-102 (2021).\n22. Blanco, E. et al. Age-associated distribution of normal B-cell and plasma cell subsets in\nperipheral blood. J Allergy Clin Immunol 141, 2208-2219 e2216 (2018).\n23. Kaetzel, C.S. The polymeric immunoglobulin receptor: bridging innate and adaptive\nimmune responses at mucosal surfaces. Immunol Rev 206, 83-99 (2005).\n24. Laumont, C.M., Banville, A.C., Gilardi, M., Hollern, D.P. & Nelson, B.H. Tumour-\ninfiltrating B cells: immunological mechanisms, clinical impact and therapeutic\nopportunities. Nat Rev Cancer 22, 414-430 (2022).\n25. Oh, Y.K., Jang, E., Paik, D.J. & Youn, J. Early Growth Response-1 Plays a Non-redundant\nRole in the Differentiation of B Cells into Plasma Cells. Immune Netw 15, 161-166 (2015).\n26. Nakamura, S. et al. Single-cell transcriptome analysis reveals functional changes in\ntumour-infiltrating B lymphocytes after chemotherapy in oesophageal squamous cell\ncarcinoma. Clin Transl Med 13, e1181 (2023).\n27. Young, M.D. & Behjati, S. SoupX removes ambient RNA contamination from droplet-\nbased single-cell RNA sequencing data. Gigascience 9 (2020).\n28. Zhang, Y., Zhang, J., Chen, X. & Yang, Z. Polymeric immunoglobulin receptor (PIGR)\nexerts oncogenic functions via activating ribosome pathway in hepatocellular carcinoma.\nInt J Med Sci 18, 364-371 (2021).\n29. Tey, S.K. et al. Patient pIgR-enriched extracellular vesicles drive cancer stemness,\ntumorigenesis and metastasis in hepatocellular carcinoma. J Hepatol 76, 883-895 (2022).\n30. Biswas, S. et al. IgA transcytosis and antigen recognition govern ovarian cancer immunity.\nNature 591, 464-470 (2021).\n31. Ocak, S. et al. Loss of polymeric immunoglobulin receptor expression is associated with\nlung tumourigenesis. Eur Respir J 39, 1171-1180 (2012).\n61"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Q5 was not appropriately addressed, as scRNA seq is feasible with FFPE as well for proper validation.\nq6 was not appropriately addressed the validation cohort still only has p16 IHc testing.\nq7 was not appropriately addressed the current calling methodology still missing the appropriate TLS markers.\nq10 not appropriately addressed\nq12 the authors highlight major error in a simple classification based on basic TNM data, that poorly reflects on the rigor in the scientific procedures followed in the team's research conduct.\nq14 is a major issue with 2 very small cohort of 4 and 6 pts using only penile primary tumors with different methodologies raises major concerns that were not appropriately addressed.\nq17 not appropriately answered as flow can be performed on all of the 90 cases no appropriate justification of the selection of only 20 samples.\nq21 also highlights another selection bias risk without a clear processing and selection of samples defined at the time of start of protocol based on the availability of tissue it seems that the authors selected a very small different cohort from the large validation cohort diluting and raising question of validity of the presented data.", "reviewer_2": "My previous questions have been nearly addressed. However, there are still some issues that require further discussion.\n1. Figures 6A-E demonstrate that the presence of IgA in the B cell non-clearance group has a pro-tumor effect in the absence of PIGR, whereas this phenomenon is not observed in IgA-/- mice. I believe this warrants discussion by the authors.\n2. The level of IgA in F4J seems to be much less than in F3J, please explain.\nMinor comments\n1. In line 298, it should be S4E, not 4E.\n2. In F4J, please clearly indicate which is HPV- and which is HPV+.", "reviewer_3": "The authors of the manuscript \"A B cell-IgA-epithelial axis enhances antitumor immunity and improves outcome in HPV-associated penile squamous cell carcinoma\" have addressed all my comments and I recommend the publication of the manuscript.", "reviewer_4": "The authors have addressed my fundamental concerns.\nEverything can get better but I believe that the human data shown in the manuscript are relevant and new enough, given the rarity of these tumors.\nThe identification of TLS formation and dominant IgA responses in penile cancers is novel and significant.\nUnfortunately, only 2 of the 3 references that the authors claim to have added are shown in the revised mansucript (response to the first paragraph). No need to change anything at this stage, but a little more rigor would have been appreciated.\nAlthough mouse data are less relevant, the improvement in the presentation of human results, and the new concepts on the immunobiology of this disease, make this study novel and exciting."}, "author_response": "REVIEWER COMMENTS (NCOMMS-24-16489B)\nReviewer #1 (Remarks to the Author):\nI appreciate the authors' detailed responses, however, several significant methodological\nconcerns remain not fully resolved:\nQ1：scRNA-seq was performed on only 13 patients, despite the existence of FFPE-based\nprotocols that, while technically limited, could have allowed expansion to a larger number of\nsamples. While I understand the technical limitations of FFPE-derived scRNA-seq, even low-\nyield data from 40-50 cases using optimized FFPE protocols could serve as partial validation,\nespecially for reproducibility of key cellular findings. I would recommend either completing the\nabove or not using the data generated from the scRNAseq for definitive conclusions.\nQ2: HPV testing was limited to p16 IHC in the majority of the cohort, despite the feasibility of\nPCR on FFPE tissue. The claim that PCR was only possible on the 30 most recent fresh\nsamples is not sufficiently justified, as FFPE blocks are commonly used for DNA-based\nPCR assays. The decision not to attempt PCR across all samples leaves room for unaddressed\nHPV heterogeneity and undermines classification fidelity.\nQ3: Risk of selection bias: The selection of subsets of patients for different assays, without a\nclear pre-defined stratification plan, introduces the risk of unconscious selection bias. For\ninstance, the 13 patients used for scRNA-seq and the 20 used for flow cytometry were drawn\nfrom different subsets, but the rationale for their inclusion remains unclear. This fragmented\ndesign challenges the generalizability of findings. A clearer explanation of how these sub-\ncohorts were selected—beyond availability—would help reduce concern about sampling bias.\nIf tissue availability was the only determinant, this should be explicitly stated with detailed\nsupplementary documentation of inclusion/exclusion criteria.\nFlow cytometry was conducted on only 20 patients, raising concerns about the\nrepresentativeness of the immune profiling. Again, given the importance of immune\ncharacterization to your study conclusions.\nQ4：Sample size and methodological inconsistency: The extremely small discovery cohorts\n(4 and 6 patients), each profiled using distinct methodologies, remain a significant weakness.\nAlthough the authors have added 3 new cases, this does not sufficiently mitigate the concerns\ngiven the continued methodological disparity between the discovery and validation cohorts.\nThe limited sample size also reduces the statistical power of exploratory analyses.\nThank you for correcting the TNM classification data. I acknowledge the effort made in\naddressing this issue.\nGiven these unresolved issues, the manuscript would significantly benefit from additional\nclarification, methodological standardization, and expanded analyses from potential other\nsites.\nWe appreciate Reviewer #1’s concerns and recognize the importance of addressing these\npoints to improve the manuscript. In response, we have provided additional explanations and\nclarifications to better justify our approach and conclusions. We hope that the rationale behind\nthese decisions will address the remaining concerns. We remain open to further suggestions\nand are happy to discuss alternative approaches if needed. Below, we present our point-by-\npoint responses to the reviewer #1' comments.\nPoint-to-point responses:\n1. Sample size and rare cancer standards\nOur cohort of 90 PSCC cases (60 archival FFPE and 30 fresh/FFPE samples) represents over\na decade of continuous collection for a rare cancer with an annual incidence of <1 per 100,000.\nThis sample size significantly exceeds the median cohort sizes reported in previously\npublished PSCC studies (see Supporting Table 1), and reflects substantial effort in the context\nof rare tumor research. Furthermore, our literature review indicates that scRNA-seq studies\nin rare cancers typically include 5–20 samples (Supporting Table 2). Our dataset of 13\nhigh-quality fresh PSCC scRNA-seq profiles—recently expanded to 15—falls well within this\nestablished range and is consistent with field standards and peer-reviewed publications.\nWe respectfully note that the expectation of performing scRNA-seq on 40–50 archival FFPE\nsamples is neither consistent with standard practice in rare cancer research nor\ntechnically feasible. The RNA integrity of our archival FFPE samples is low (median DV200\n< 20%), well below the generally accepted threshold (DV200 > 30%) required for successful\nlibrary construction. Under such conditions, attempting high-throughput transcriptomic assays\nwould likely yield poor-quality data and constitute an inefficient and potentially wasteful\nuse of scarce archival tissue—raising both scientific and ethical concerns.\n2\n2. HPV detection strategy is evidence-based and consistent with field standards\nOur classification of HPV status adheres to the guidelines of the International Society of\nUrological Pathology (ISUP) (PMID: 29105175; PMID: 32235153), which recommend p16\nIHC as a validated surrogate for HPV positivity in PSCC. This recommendation is supported\nby extensive validation studies (e.g. lmmunity 2020 PMID: 31924475; Nature Communications\n2021, PMlD: 34099645). To further validate this approach, we performed confirmatory HPV\nPCR on 30 fresh tumor samples. This validation analysis has since been expanded to 35\ncases, including the five most recent samples (P91–P95), all of which showed complete\nconcordance with p16 IHC status (Supporting Fig.1). These findings reinforce the reliability\nof p16 IHC classification strategy.\nWhile Reviewer #1 suggests that PCR should have been conducted across all archival FFPE\nsamples, we concerned this expectation overlooks well-documented technical limitations.\nHistorical FFPE blocks are often compromised by DNA degradation and low yield, making\nthem unsuitable for sensitive PCR-based detection. According to College of American\nPathologists (CAP) guidance and multiple prior studies, HPV PCR is not recommended when\nDNA integrity is insufficient, as this may result in a high false-negative rate (PMID:37351572,\nPMID:28812131, PMID:26514313, PMID:25421801). Applying PCR under such conditions\nwould introduce systematic misclassification and reduce overall accuracy (Supporting\nFig.2).\nDespite these significant concerns, and in an effort to accommodate the reviewer’s\nexpectations, we still proceeded with HPV PCR on archival FFPE samples. As anticipated,\nwe observed 5 instances of p16-positive but PCR-negative cases in 60 archival FFPE\nsamples (Supporting Fig.3). Given the known limitations of PCR sensitivity in degraded FFPE\ntissue, we believe these discrepancies most likely represent false negatives due to technical\nartifacts, not inaccuracies in p16-based classification. Nonetheless, if the reviewer prefers a\nmore conservative analytical approach, we are willing to exclude the five discordant FFPE\nsamples (p16-positive but PCR-negative) from the analysis. This adjustment would reduce\nthe total cohort size from 90 to 85 cases. We have reanalyzed the data accordingly, and\nimportantly, the key statistical conclusions remain unchanged (revised Figures in\nmanuscript)\nHPV DNA positivity (~180 bp band) showed complete concordance with p16 IHC status across\nall cases.\nFigure 2. HPV PCR shows false-negative results in archival FFPE samples compared to\np16 IHC and HPV E6/E7 staining\n(A) Representative images of HPV status determined by p16 IHC and HPV16/18 E6/E7\nimmunofluorescence. P34 and P41 show strong p16 expression and positive E6/E7 signals,\nconsistent with HPV-positive status. Ctrl: HPV-negative case with no p16 or E6/E7 staining.\n(B) Sequences of HPV-specific universal primers SPF1/GP6++ used for PCR detection.\n(C) HPV DNA detection by PCR in the same four cases. While the positive control shows a\n~180 bp band indicating HPV positivity, no bands are detected in P34 and P41, despite their\nconfirmed HPV-positive status by p16 and E6/E7 staining—highlighting false-negative PCR\nresults likely due to DNA degradation in FFPE samples.\nHPV DNA PCR was performed on 60 archival FFPE samples to evaluate concordance with\np16 IHC. While the majority of samples showed concordant results, five cases were p16-\npositive but PCR-negative, likely reflecting false-negative PCR results due to degraded DNA\nquality inherent in FFPE tissue.\n3. No selection bias: transparent, time-sequential cohort accrual and assay allocation\nWe respectfully clarify that our study design did not involve any form of pre-selection\nbased on biological characteristics or expected outcomes. All fresh PSCC samples were\ncollected prospectively and consecutively from 2021 to 2024 and have since been\nexpanded to 35 cases, including the most recent additions (P91-P95). Tissue allocation for\ndownstream assays was determined solely by real-time sample availability and the\ntechnical compatibility of the assays (Supporting Fig. 4)\nFigure 4. Time-sequential accrual and assay allocation workflow of the PSCC cohort.\nconsumed for scRNA-seq, which requires high tissue input and leaves no residual material.\nAs a result, these samples could not be used for subsequent assays such as flow cytometry.\nAs additional cases were collected over the following years, each was directly processed for\nflow cytometry, yielding 20 total cases. All were analyzed in the order received. This\nsequential and availability-driven workflow is standard practice in rare cancer research,\nwhere tissue is scarce and assays may be mutually exclusive.\nWhile the scRNA-seq and flow cytometry cohorts do not overlap, we respectfully disagree with\nthe characterization of our study as a “fragmented design” that undermines generalizability.\nOn the contrary, these two datasets represent independent, sequentially accrued subsets\nfrom the same clinical population, analyzed using complementary technologies. This structure\nprovides a form of natural replication across time and platform, and the fact that both\nindependently identify TIL-B-cell enrichment in HPV+ tumors strongly supports the\nreproducibility of this finding.\nImportantly, this observation was further validated in a larger cohort (n=85, excluding 5\ndiscordant p16+/PCR− cases) using immunohistochemistry, confirming the same immune\nsignature. Thus, our conclusion is not only reproducible but triangulated across three\northogonal, methodologically distinct platforms—scRNA-seq, flow cytometry, and IHC—\neach reinforcing the others. This multi-platform concordance provides rigorous cross-\nvalidation, and directly addresses any concern about individual cohort size or assay-specific\nlimitations.\nIn this context, we would like to clarify that the flow cytometry data were never presented as\nthe sole or primary basis for our conclusions. Rather, they served as a supportive\ncomponent consistent with findings from scRNA-seq and IHC. If the reviewer continues to\nregard the sample size of the flow cytometry cohort as insufficient, we remain open to\nremoving this dataset from the revised manuscript, if absolutely necessary. However, we\nemphasize that doing so would not alter the key conclusion of this study—namely, the\nenrichment of TIL-B cells in HPV+ tumors—which is robustly supported by multiple\nindependent and orthogonal approaches.\n4. Use of 3′ and 5′ scRNA-seq platforms does not constitute methodological\nheterogeneity\nWe respectfully clarify that the use of both 3′ and 5′ single-cell RNA-seq chemistries in our\nstudy is methodologically sound and widely accepted in the field. While the two platforms\ndiffer in transcript capture location, both reliably generate gene-level expression profiles\nsuitable for cell-type annotation, clustering, and differential expression analysis (PMID:\n39675380; PMID: 33472597). Integration of 3′ and 5′ datasets has become a standard,\nvalidated practice in single-cell transcriptomics, and has been successfully applied in\nnumerous high-impact studies (PMID: 33657410; PMID: 32079746), using well-established\ntools such as Seurat v3+, Harmony, and LIGER.\nrecently expanded our PSCC cohort from 13 to 15 high-quality scRNA-seq profiles. As an\nalternative analysis, we also evaluated a subset of 11 cases generated exclusively using\n3′ scRNA-seq, excluding the 4 cases based on 5′ scRNA-seq. Importantly, this 3′-only\nanalysis fully recapitulated our original findings, including statistically robust differences\nbetween HPV+ and HPV− groups (Supporting Fig.5). These results further support the\nreproducibility and robustness of our conclusions, regardless of sequencing platform. We\nare happy to present the 3′-only analysis in the revised manuscript, should the reviewer or\nEditor prefer this presentation.\nFigure 5. Reproducibility of statistically significant enrichment of TIL-Bs in HPV+ PSCC\nusing 3′ scRNA-seq–only analysis\n(A) UMAP plot of 34,119 immune cells from 3′ scRNA-seq profiles (HPV−, n = 5; HPV+, n =\n6), visualized across 12 transcriptionally distinct clusters and color-coded by cell type.\n(B) Relative abundance of major immune cell subsets (TIL-Bs, TIL-Ts, NK cells, and myeloid\ncells) across individual HPV− and HPV+ PSCC samples.\n(C) Bar graph comparing the group-level proportions of immune cell types between HPV− and\nHPV+ tumors.\nThis 3′-only analysis—excluding four 5′ samples—recapitulates the original results and\ndemonstrates statistically significant enrichment of TIL-Bs in HPV+ PSCC.\nStatistical analysis was performed using a two-sided unpaired Student’s t-test. Data are\npresented as mean ± s.e.m. with significance indicated as: ns, not significant; **P < 0.01.\nMy previous questions have been nearly addressed. However, there is still one issue that\nrequire further discussion.\nQuestion: The apparent discrepancy between the representative clonogenic assay images in\nFigure 5D and the corresponding quantitative analysis warrants further clarification. The\nPIGR-overexpression group appears to demonstrate enhanced colony formation in the\ngraphical presentation. To support the conclusion of PIGR-overexpressing PSCC primary\ncells showed reduced colony formation, it is recommended that the authors carefully verify the\nexperimental replicates and provide all images in statistics for review.\nResponse:\nWe thank the reviewer for pointing out the discrepancy between the representative images\nand the quantitative data in Figure 5D. Upon careful review, we found that the issue was\nprimarily due to a reduction in image resolution, which inadvertently gave a misleading\nimpression of colony density in the PIGR-overexpression group. After correcting this and\npresenting the image at higher resolution, the visual presentation is now consistent with the\nquantified results. As shown in the revised figure, PIGR-overexpressing PSCC primary cells\nexhibit significantly reduced colony formation (Figure below). To ensure full transparency, we\nhave also provided the original, unprocessed images used for quantification in the\nsupplementary materials for the reviewer’s reference.\nFigure 6. The initial low-resolution image inadvertently gave a misleading visual impression\nof colony density in the PIGR-overexpression group (Left). Upon reprocessing at higher\nresolution, the image more accurately reflects the reduced colony formation phenotype (Right).\nThis correction addresses the artifact introduced by image compression and improves the\nconsistency between visual and quantitative data.\nPIGR-overexpressing PSCC primary cells exhibit significantly reduced colony formation.\n9"}, {"round": 3, "reviewer_comments": {"reviewer_6": "I appreciate the authors' and reviewers’ tremendous efforts and do not have further comments for the authors to address.\nThank you for inviting me as the arbitrating reviewer."}, "author_response": ""}]

Decision letter