Clinical information and Sequencing materials of 119 patients diagnosed with esophageal squamous cell carcinoma (ESCC) were collected from Anyang Cancer hospital as described previously (33). Briefly, similar sized tumor and adjacent normal samples from each patient were used to extract RNA using Invitrogen’s TRIzol Regents according to the manufacturer’s instructions. After quantification with Agilent 2100 Bioanalyzer (Agilent RNA 6000 Nano Kit), 1 μg RNA was used to construct the sequencing library following the method provided by VAHTS^® Total RNA-seq (H/M/R) Library Prep Kit for Illumina^®. Quantified libraries were sequenced using Illumina X Ten platform (BGI) with paired-end 150bp read length and at least 10 million paired-end 150bp reads for either tumor or matched adjacent normal tissue samples was generated. The raw sequencing data has been deposited at the National Genomics Data Center of China (https://bigd.big.ac.cn/) (Bioproject Access ID: PRJCA001577).

The gene expression values were calculated and represented as log2(FPKM +1). Briefly, raw sequencing data were first quantified against the indexed GRCh37 genome downloaded from Ensembl (34) using Salmon (version 0.9.0) (35). The quantified gene counts were loaded in an R 3.5 environment and genes without expression value in 50% of the samples were dropped. After that, the gene counts were normalized using the CQN (Conditional Quantile Normalization) pipeline (36) and the FPKM value was calculated using the normalized gene counts.

MixCR Software (20) was used to extract CDR3 repertoires from Fastq files. To minimize the technical error and maximize the confidence of comparative analysis, only samples with at least 500 reads mapping the Ig regions and CDR3s with clone counts ≥ 3 were used for further analysis. For V gene and J gene usage analysis, only one clone with an unambiguously mapped V or J gene was taken into account. The total clone abundance (N) was calculated as the sum of the clone counts of each clone. The clone fraction of each V gene was calculated using the clone counts of the V gene divided by N. To calculate rearrangements, we only considered the V J subgroups such as IGVH1, IGVH3, etc. The proportion of one V-J rearrangement was calculated as the sum of the clone fraction of each clone. When isotyping, clones mapped to IGHG1, IGHG2, IGHG3, and IGHG4 were counted as IGHG, while clones mapped to IGHA1 and IGHA2 were recognized as IGHA. We also classified each subclass of IGHG and IGHA, counting those clones that only mapped to one of the subclasses.

To quantify the diversity of the BCR repertoires, clonality was calculated as 1 – normalized Shannon-Wiener index. A value close to 0 indicates an even distribution of clone fraction and more diverse BCRs. On the contrary, a value close to 1 indicates BCRs are more oligoclonal.

According to the gene expression profiles of different lymphocytes, CIBERSORT (37) can deconvolve 22 immune cell subtypes. The detailed gene list used to define each immune cell type can be found at https://cibersort.stanford.edu. Firstly, we built a gene expression matrix as the CIBERSORT input mixture file. Each row of the matrix represents the expression (FPKM values) of a given gene and each data column consists of the expression profile for a single sample. Secondly, we updated the input mixture file to https://cibersort.stanford.edu/. Finally, the bioinformatic analysis was conducted online with the LM22 signature gene file (500 permutations).

Survival analysis was undertaken using the survival (https://github.com/therneau/survival) and survminer (https://github.com/kassambara/survminer) packages in the R 3.5 environment. The cut-off value to split patients into two groups was determined using the surv_cutpoint function from the survminer package. All survival plots were generated using the Kaplan-Meier estimator. Plots were created with ggsurvplot function. Multivariable analysis was performed with Cox proportional hazard regression.

As tumor infiltrating immune cells constitute only a small part of the tumor tissues, sequencing is key for discovery of TCR or BCR genes. In our dataset, at least 50 million paired-end reads were collected for each tumor or normal adjacent sample ( Figure 1A ). The number of reads successfully mapped to TCR or BCR regions ranged from 720 to 1.1 million without bias between tumor and related adjacent normal samples ( Figure 1B ). Notably, 97.84% ( ± 3.6) of aligned reads mapped to BCRs, and only a handful of reads mapped to TCRs ( Supplementary Table S1 and Figure 1C ). Previous single cell RNA-seq investigating single cell TCR profiling (15) revealed that TCR clone type compositions and proliferative cell proportions were highly diverse across different T cell subtypes. For example,CD4+ T cell subtypes had a low degree of clonal sharing and clone size. The diversity of TCR and BCR identified may reflect the differential immune cell infiltration in individual patient TMEs. In order to obtain global profiles of B cell and T cell-mediated immune responses in ESCC tumor, we further deconvolved the tumor-infiltrating immune cells using quantified gene expression data with CIBERSORT (37). ESCC was characterized by a predominant plasma cell and resting memory CD4 T cell signature, but low infiltration of activated T cells ( Figure 1D ). Therefore, the low clone counts of TCRs in an immune suppressive ESCC tumor environment may be the cause of biased alignment of TCR and BCR reads. The insufficient T cell immunity and unique signature of B cell infiltrations in ESCC promoted us to further investigate the role of B cell immunity in disease progression and control.

To minimize the effect of clones generated from technical error introduced by sequencing, we only took into account those CDR3 (third complementarity determining region of the heavy chain) clone records with no ambiguous mapping and at least 3 supporting reads for V J usage analysis. A total of 71 or 70 IGHV genes were identified in tumor and normal samples, respectively. The IGHV gene usage is very biased as the top five genes account for about 37.14% in tumor samples: IGHV3-23 (10.71%), IGHV1-18 (7.73%), IGHV4-39 (7.64%), IGHV3-21 (6.14%) and IGHV5-51 (4.92%) ( Figure 2A ). When compared to the normal samples, we observed a higher proportion of IGHV3-21 (median usage: 6.14% vs 4.09% (tumor vs normal tissue) student t-test p: 0.00019) and lower proportion of IGHV3-74 usage in tumor samples (median usage: 2.23% vs 5.49%, (tumor vs normal tissue) student t-test p: 1.02e-8) ( Figures 2B, C ). IGHJ4 and IGHJ6 were the most used J segments, and the usage of 6 J segment genes was consistent between tumor and normal tissues ( Figure 2D ). Finally, we investigated the usage of V genes in IGK and IGL chains. Although similar V gene usage distribution between tumor and normal tissue were observed, the usage of IGK and IGL V genes was also greatly biased. There were 31 IGK V genes and 46 IGL V genes in the ESCC tumor samples and the usage proportion of the top 5 V genes was 59.76% and 48.48%, respectively ( Supplementary Figure S1 ).

Because of the limited number of sequencing reads mapped to the BCRs compared to the diversity of the rearrangement, we only studied V gene subgroups to increase the confidence of analysis. IGHV3 (tumor: 40.3%, normal: 44.4%) and IGHV4 (tumor: 23.1%, normal: 23.6%) were the most frequently used heavy chain V genes ( Figure 3A ). IGHV3/IGHJ4 rearrangements are most common in ESCC tumors yet showed significantly lower usage compared to normal samples with a median frequency of 21.5% vs 24.6% (tumor vs normal tissue; student t-test p: 3.76e-06, Figure 3B ).

In the 119 tumor samples, 105614 unique IGH CDR3 amino acids were found. Neither the CDR3 length distribution nor the clone abundance was different between tumor and normal samples, yet wide spectrums of CDR3 clone counts were observed among individuals ( Supplementary Figure S2 , Figures 3C, D ). Only 4.5% of the IGH CDR3 clones were shared between tumor and normal samples ( Figure 3D ), which is much lower than that of IGL (14.5%) and IGK (35.8%) chains. We further calculated Jaccard Index for each paired tumor and normal sample. This index ranges from 0 to 1 and the higher the percentage, the more similar the two populations. Results showed that the B cells in tumors and the matched normal samples were almost two distinct populations with extremely low similarities as the median IGH Jaccard index was 0.006. In line with the amount of shared CDR3 clones, the Jaccard index for IGK and IGL chains was significantly higher than the heavy chain (Student t-test p< 0.0001, Supplementary Figure S2 ).

We calculated the 1 - normalized Shannon-wiener index of all IGH Chain CDR3s considering both CDR3 contents and frequency as a representation of clonality. Intriguingly, tumor samples showed a more oligoclonal profile of B cells compared to matched adjacent normal samples (Student t-test p: 5.41e-18) as the high-frequency clonotypes were exclusively observed in tumor samples ( Figure 4A ). When one base mismatch or deletion was defined as the main type of CDR3 somatic hypermutation (SHM), the median mutation ratio for tumor samples is 4.7%, significantly higher than that of normal samples (Wilcoxon ranksum p: 0.01). There was no association of SHM status between tumors and the corresponding normal samples ( Supplementary Figure S3 ). In agreement with findings in normal samples and a previous report (23), the mutational load in tumor samples was positively associated with clone abundance (Spearman rho: 0.956) and the majority of mutations were synonymous ( Figure 4B ).

We defined CDR3 clones mapped to IGHG1, IGHG2, IGHG3 and IGHG4 as IGHG repertoires, IGHA1 and IGHA2 as IGHA repertories. To minimize the measurement deviation caused by sequencing depth, clonotypes per one million sequencing reads (CPM) (38) was introduced to quantify each Ig subclass. CPMs were highly consistent with RNA-seq-quantified FKPMs for all Ig subclasses ( Supplementary Figure S4 ). Compared with normal samples, TIBs in ESCC showed a clear predominance of IGHG in most cases. Nonetheless, the CPM distribution for each IgG subclass showed differences and IgG1 quantification was much higher than other IgGs ( Figures 4C, D ). Aside from the constraints of RNA-seq data and expression differences, this may also be caused by subsequent class switch recombination (CSR) (23). When only clones mapped to two different subclasses were taken into account, we found that 18322 IGHG4 clones co-existed with other subclasses, and IgG2/4 (7.76%), IgG3/1 (3.14%) ranked as the top two combinations ( Figure 4E ).

We first examined the possible prognostic value of the expression level of B cell antibodies. Patients were split into low and high expression groups according to the gene expression distribution, and at least 20 patients were included in each group. Kaplan-Meier analysis showed that higher expression of IGH (a sum of the FPKM values of IGHG1, IGHG2, IGHG3, IGHAG4, IGHA1, IGHA2, IGHD, and IGHM genes) is significantly associated with longer overall survival with a hazard ratio (HR) of 0.499 (log-rank p: 0.014, Figure 5A ). The expression level of IGHG2, IGHG3, and IGHG (a sum of the IGHG1, IGHG2, IGHG3, IGHAG4) also had a positive effect on survival ( Supplementary Figure S5 and Figure 5B ).

Secondly, we examined the prognostic potential of B cell clonality and class-switch recombination (CSR). Kaplan-Meier survival curve analysis showed that higher clonality (lower diversity) was negatively associated with prolonged survival with a HR of 1.97 (log-rank p: 0.012, Figure 5A ). As to CSR, we also confirmed that high levels of CSR are significantly associated with longer survival (log-rank p: 0.002, Figure 5A ). Regardless of the important role of IgG1/IgG3 isotypes in initiating tumor-specific antibody-dependent cellular cytotoxicity (ADCC), high IGHG3/1 class switch events were a neutral survival parameter in ESCC. On the other hand, IGHG2/4 and IGHG1/2 class switch events are significantly associated with better prognosis (log-rank p: 0.004, 0.034) ( Supplementary Figure S6 ). A multivariate analysis using Cox proportional-hazard regression with adjustment for gender, age, smoking, drinking, tumor stage and infiltration was further conducted and confirmed the prognostic association of the above B cell characteristics ( Figure 5B ).

The main TIBs resolved by CIBERSORT were plasma cells. We further assessed the prognostic value of the plasma cell marker CD138 and non-plasma B cell marker MS4A1 (encoding CD20). In line with another study (39), high expression of CD138 was strongly beneficial for ESCC patient survival (log-rank p: 0.003). This association became more significant when we examined the ratio of IGH to MS4A1, which reflects the relative abundance of antibody-producing plasma cells compared to CD20-positive non-plasma B cells (log-rank p: 0.001, Figure 6A and Figure 5B ). Since the Ig isotypes showed varied prognostic values in ESCC, we further assessed the ratio of isotype-specific plasma cells to non-plasma cells with multivariable Cox proportional-hazards regression and found that IgG1/3/4 expressing plasma cells had a neutral effect on prognosis, but high levels of IgG2-expressing plasma cells, determined by the ratio of IGHG2 to MS4A1(IGHG2 ratio), was strongly associated with better survival (HR: 0.40, Figure 6B ). We further explored the association between IGHG2 ratio and pathological characteristics and found that the IGHG2 ratio significantly decreased in III/IV stage patients compared with stage I or II patients ( Figure 6C ). High expression of IGHG2 is also a positive survival parameter in different solid tumors when validated using TCGA (The Cancer Genome Atlas) data, including head and neck squamous cell carcinomas (TCGA-HNSCC) and breast invasive carcinoma (TCGA-BRCA) ( Supplementary Figure S7 ). Taken together, this data suggests that IgG2-expressing plasma cells play an important role in B cell anti-tumor immunity in ESCC.

According to differences in gene expression profiles, we recently classified ESCC into four molecular subtypes: differentiated, characterized by keratinocyte differentiation and epidermis development; immunogenic, characterized by B-cell surface markers and T-cell chemokine ligands; metabolic, defined by upregulation of genes involved in drug metabolism by cytochrome P450 and retinol metabolism; and stemness, associated with stem cell markers and the Wnt signaling modulator SFRP1. Each molecular subtype was found to contain a very different TMEs and subtypes are strongly associated with prognosis (33).

In line with gene expression profiles of the molecular subtypes (33), TIBs were most abundant in the immunogenic subtype (median CDR3 clones: 65351), whereas the stemness subtype (median: 4888) was much less B cell infiltrated (Kruskal–Wallis test p: 1.2e-14). We also observed the IGHV4-39 usage proportion (median: 3.63%) was significantly decreased in stemness samples (Kruskal–Wallis test p: 0.028). No statistical differences were found for IGHV3-74 frequency, clonality, IgG2, or IgA CSR ratio among the four molecular subtypes, yet the stemness samples had a higher IgG3/1 ratio compared to other subtypes (Kruskal–Wallis test p: 0.0018) ( Figure 7A ). In addition, the immunogenic samples had the highest expression of MS4A1 but lowest CD138 to MS4A1 ratio, suggesting that the TIBs in the immunogenic subtype are mainly CD20 positive non-plasma cells. Although the differentiated and metabolic subtypes of ESCC have higher CD138 expression than immunogenic and stemness subtypes, the CD138 to MS4A1 ratio of stemness was comparable with differentiated and metabolic subtypes ( Figure 7B ).

When we investigated the IGHG2 to MS4A1 ratio between subtypes, no significant difference was observed. However, multivariable Cox analysis confirmed that IgG2 level could further stratify the high-risk patients within immunogenic (HR: 0.319, log-rank p: 0.03), differentiated (HR: 0.142, log-rank p:0.021) and metabolic subtypes (HR: 0.057, log-rank p: 0.00013) ( Figure 7C ). As to the stemness subtype, which has the worst clinical outcome, the IgG2 plasma cell level showed no association with better survival (log-rank p: 0.3). Another multivariable Cox analysis was further conducted for non-stemness patients, resulting in a hazard ratio of 0.253 (log-rank p: 0.0006, Figure 7C ).

Thus, combining the molecular subtyping and IgG2 expressing plasma cells in the TME could greatly benefit ESCC patient screening by determining their relative risk.