Fresh tumor, adjacent non-tumor lung tissue, and PBMC were obtained from early-stage LUAD patients (n = 49) undergoing surgical resection at the James Graham Brown Cancer Center, University of Louisville. Patients were treatment naive or, in the case of three patients, had not received chemo/radiation treatment at least 65 days prior to surgery. Samples were anonymously coded in accordance with the Declaration of Helsinki; written informed consent was obtained under protocols approved by the University of Louisville School of Medicine IRB. Cells were isolated, counted, resuspended in FBS/10% DMSO, and stored at –150°C.

Fresh tumor and adjacent tissues (n = 12) were minced and digested in RPMI 1640 with 2% FBS, type IV collagenase (1 µg/ml), and hyaluronidase (10 ng/ml) for 20–40 min at 37°C. Single-cell suspensions were stained with fluorescently labeled mAbs (CD45, CD19, CD69, CD27, CD38, viability dye) and analyzed on a BD Canto II flow cytometer.

Human lung cancer and distant tissues were snap-frozen in OCT and stored at –80°C. Eight-µm sections were fixed with cold acetone for 15 min and air-dried for 30 min. Slides were blocked with 20% FBS in PBS for 1 h, stained overnight at 4°C with anti-CD20-PE (1:100), CD3-FITC (1:100), and DAPI, and imaged on a Nikon A1R confocal microscope.

Preserved cells were thawed quickly in a 37°C water bath and washed three times in pre-warmed complete RPMI 1640 media and then counted to determine proper input. Cell aliquots were then resuspended in serum free RPMI and incubated with Cisplatin (Fluidigm; Cat#201064, 1:1000) for 5 minutes at room temperature. After incubation, cells were washed with 5X volume of pre-warmed complete RPMI 1640, then re-suspended and washed in 4mL Maxpar Cell Staining Buffer (Fluidigm, Cat#201068). Cells were incubated for 10 minutes at room temperature with 2uL of Human TruStain FcX (Biolegend) and stained with a panel of 44 CyTOF antibodies for 30 minutes. CyTOF antibodies were either purchased pre-conjugated (Fluidigm) or conjugated in house using either Maxpar MCP9 antibody labeling kits or Maxpar X8 antibody labeling kits (Supplementary Table 1). Cells were washed with Maxper Cell Staining Buffer and fixed with 1.6% Formaldehyde (Thermofisher, Cat#28906) in Maxpar PBS (Fluidigm, Cat#201058) for 10 minutes at room temperature. Then cells were re-suspended with 1mL of Cell-ID™ Intercalator-Ir—125 µM at 1:8000 in Maxpar Fix and Perm Buffer (Fluidigm, Cat#201067) and incubated at 4°C. The following morning, cells were washed with 2mL Maxpar Cell Staining Buffer once and 1mL of Maxpar Cell Acquisition solution (Fluidigm, Cat#201241) once. Pelleted cells were kept on ice until analysis. Cells were then resuspended in EQ Four Element Calibration Beads (Fluidigm, Cat#201078)/Maxpar Cell Acquisition solution at 1:10 ratio with cell concentration around 1.0x10⁶ cells/mL. Cells were filtered before acquisition on a Helios (Fluidigm) at approximately 500 events/second. Standard gating analysis identifying the number of live immune cells (CD45+; mean=82,226 cells per patient), T cells (CD3+CD4+, CD3+CD8+; mean=44,616 cells per patient) and B cells (CD19+; mean=9,344 cells per patient) was conducted through FlowJo software (BD Life Sciences) for subsequent clustering and differential analysis.

CyTOF data was analyzed using the CyTOF Workflow (25, 26). Briefly, FlowJo workspace files were imported and parsed using functions from flowWorkspace and CytoML, then stored as a single cell experiment object after arcsinh transformation through CATALYST (25, 27). Unsupervised cell clustering was conducted using FlowSOM and ConsensusClusterPlus in two parts (25, 26). Initial clustering was done by first compiling all events from each sample and using primary lineage markers CD19, CD3, CD4, CD8, and CD11b to generate 20 unique clusters. Clusters were then labeled as CD19+ B cells, CD3+CD4+ T cells, CD3+CD8+ T cells, or CD11b+ myeloid cells and merged for use in secondary clustering. CD56+CD16+ clusters or CD3+TCRgd+ clusters were not merged and labeled NK or TCRgd T cells, respectively. Each merged cluster was then further dissected into 15 total subsets through subclustering based on median expression of both lineage and selected functional markers (Table 1).

RNA was extracted from matched PBMC, tumor, and normal tissue samples of 29 patients from the CyTOF cohort, using RNeasy Mini Kit (Qiagen) according to manufacturer instructions. Yields were measured using Qubit RNA High Sensitivity Kit (ThermoFisher Scientific) and integrity was measured using the Bioanalyzer RNA Nano Kit (Agilent Technologies) and stored at -80C. RNA was prepared for sequencing using the SMARTer Human BCR IgG IgM H/K/L Profiling Kit (Takara Bio USA) according to manufacturer instructions with no modifications. Briefly, approximately 100ng of RNA was used for PBMC samples and 400ng for tumor and adjacent tissues. IgG and IgM primers were used to amplify heavy chain transcripts through a 5’RACE reaction. Library quality was measured using the 2100 Bioanalyzer High Sensitivity DNA Assay Kit (Agilent Technologies) and the Qubit 3.0 Fluorometer dsDNA High Sensitivity Assay Kit. Sequencing was done using the 600-cycle MiSeq Reagent Kit v3 (Illumina) on the MiSeq platform using 300bp paired-end reads. Samples were equimolar pooled into multiplexes of 24, ensuring all tissues and isotypes from each patient were grouped together, when possible. FASTQ reads were generated using the associated DRAGEN software package (Illumina).

Targeted amplification of IgG heavy chain transcripts was done using a novel near-full length sequencing method as described previously (28). Briefly, extracted RNA from a subset of patient tumor samples (n=21) used in CyTOF and traditional AIRR-seq was thawed on ice. Approximately 100–600 ng of RNA was used per sample and converted to first strand complementary DNA (cDNA) using the SMARTer RACE 5’/3’ Kit (Takara Bio USA) according to manufacturer instructions. A custom template switch oligonucleotide containing a unique molecular identifier (5’ TSO-UMI) was used for template switch during cDNA synthesis. PCR amplification of heavy chain transcripts was done using 10uL 5X PrimeSTAR GXL buffer, 4uL GXL dNTP mixture, 28uL PCR-grade water, 1uL PrimeSTAR GXL Polymerase, 1uL 10X UPM from the SMARTer RACE 5’/3’ kit and 1uL of appropriate barcoded IgG-specific 3’ primers (Supplementary Table 2). Amplification conditions for full-length IgG were done as previously described (26). Amplified products were purified using a 1x1 cleanup with ProNex magnetic beads (Promega) and quantified using Qubit dsDNA HS assay (ThermoFisher Scientific) and length was evaluated using the Fragment Analyzer Genomic DNA HS assay (Agilent). Amplifications were equimolar pooled into eight-plexes and prepared for sequencing using SMRTbell Express Template Prep Kit 2.0 (Pacific Biosciences) according to the “Procedure and Checklist for Iso-Seq Express Template for Sequel and Sequel IIe systems” protocol starting at “DNA Damage Repair” with modifications as previously described. Final libraries were quantified using the Qubit dsDNA HS assay (ThermoFisher Scientific) and quality was evaluated using the Fragment Analyzer Genomic DNA HS assay (Agilent). Sequencing of each sample pool was performed on one SMRTcell 8M using primer v4 and polymerase v2.1 on the Sequel IIe system with 30hr movies. Reads were demultiplexed and high-fidelity circular consensus sequences (“HiFi reads”) were generated on the instrument for downstream analysis.

Both short-read AIRR-seq and FLAIRR-seq datasets were processed and analyzed using tools and packages from the Immcantation Portal (29). Quality control processing of AIRR-seq datasets was done using pRESTO (30) functions, “FilterSeq trimqual” to remove bases with <Q20 read quality and “FilterSeq length” to remove reads shorter than 125 bp. Primer sequences were identified with an error rate of 0.2 and noted in FASTQ headers using “MaskPrimers align.” Reads with the same 12 bp UMI were aligned using “AlignSets muscle” and a consensus sequence for each UMI was generated using “BuildConsensus.” Mate pairing was done with “AssemblePairs sequential” using reference genome-guided alignment with a minimum 5bp overlap and reads <400bp were removed using “FilterSeq.” Reads with the same UMI were collapsed using “CollapseSeq” and reads with <2 supporting sequences were removed using “SplitSeq group.” For pRESTO processing of FLAIRR-seq datasets, HiFi reads with <Q20 read quality and/or reads with <750 bp were removed using “FilterSeq quality” and “FilterSeq length.” Primers and a 22bp UMI were identified using “MaskPrimers align” with an error rate of 0.3 and reads with the same UMI were grouped and aligned using “AlignSets muscle.” FLAIRR-seq reads are done as a single molecule so mate pairing was not required. Subsequent quality control and filtration of consensus reads was done as described above. Filtered reads for both datasets were further processed using the Change-O (29) toolkit. For gene assignment, we used "IgBLAST" to align sequences against immunoblobulin heavy chain V (IGHV), D (IGHD), and J (IGHJ) alleles downloaded from the International IMmunoGeneTics Information System (IMGT) database (downloaded Febuary 21, 2022). Clonal clustering was performed with “DefineClones.py” followed by germline reconstruction and conversion using “CreateGermlines.py.” To define immunoglobulin heavy chain constant (IGHC) genes from FLAIRR-seq data, reads were filtered by IGHC length (900-1100bp) and aligned to chromsome 14 hg38 reference using minimap2 (31) and SAMtools (32) to generate sorted and indexed bam files. WhatsHap (33) was used to identify, genotype, and phase single nucleotide variants (SNVs), then reads from each gene were clustered using CD-HIT (34) with a 100% identity threshold. A single representative read from clusters representing at least 5% of the total reads per IGHC gene were aligned to an in-house IGHC reference dataset using BLAST (35) to confirm the IGHC gene assignments and infer alleles.

All repertoire-derived metrics were analyzed in R (v4.3.1) using functions from Alakazam (29), dplyr, and stats. Simpson’s diversity index was computed using “alphaDiversity” and 95% confidence intervals were obtained using 200 bootstrap replicates implemented within the function. Repertoire polarization was calculated using “countClones” to identify the minimum number of clones representing 80% of sequences, normalized to the total clone count. SHM frequencies were computed with "observedMutations" to calculate the number of abse changes between the full sequence alignment or specified region (complementarity determining regions and framework regions) and the constructed germline alignments.

IGHV, IGHD, and IGHJ gene usage frequencies were calculated using “countGenes.” To assess tumor enrichment, the frequency of each gene in PBMC IgM was subtracted from its frequency in tumor IgM or tumor IgG for each individual. Clone sizes for enriched versus depleted genes were compared by normalizing clone counts to total productive sequences per isotype and tissue. Principal component analysis (PCA) of gene-usage frequencies was conducted with “prcomp.” Differences in repertoire structure were quantified by calculating Euclidean distances along PC1–PC2 between (1): each sample and its nearest neighbor of the same tissue type, and (2) paired PBMC–tumor samples from the same patient.

Statistical comparisons between paired tissues were performed using paired t-tests, and differences between patient groups were compared using unpaired t-tests. When multiple genes, groups, or subisotypes were tested simultaneously, p-values were adjusted using Bonferroni correction to control for multiple hypothesis testing and significance thresholds after correction are reported in the text or figure legends.

All human samples were obtained from LUAD patients under protocols approved by the University of Louisville School of Medicine Institutional Review Board, in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants prior to sample collection. For all patient data and images, consent was obtained for use in research, and records of consent have been retained. No identifying facial material is included.

We analyzed tumor and matched non-tumoral adjacent lung tissue from a pilot cohort of 12 patients with LUAD (Supplementary Table 3). Flow cytometry revealed a significantly higher frequency of CD19+ B cells within tumors compared to adjacent tissue (Supplementary Figures 1A, B, paired t-test; p=8.02x10^-5). Within the B-cell compartment, tumors were enriched for CD19+CD27+CD38+ (p=1.72x10^-5) and CD19+CD69+ (p=0.033) subsets, whereas CD19+CXCR5+ (p=0.002) B cells were reduced relative to adjacent tissue (Supplementary Figure 1C; paired t-tests). Additionally, confocal microscopy of a subset of these patients revealed increased co-localization of CD3+ and CD19+ cells in tumor compared to adjacent tissue (Supplementary Figure 1D), consistent with previous studies (8, 11, 19, 36). Together, these initial findings highlighted an altered B cell population in LUAD tumors and underscored the need to investigate at higher resolution how B cell subsets and their expressed antibodies may shape the tumor microenvironment.

We collected peripheral blood, tumor tissue, and non-tumoral adjacent lung tissue (adjacent) from 48 LUAD patients (I n=28; II n=12; III n=7; IV n=1), divided equally between male and female (Table 2). All patients were treatment naive or had not received treatment at least 65 days prior to surgery. To characterize cell populations, we performed CyTOF analysis using a panel of 44 immune cell surface markers (Supplementary Table 1), allowing for resolution of major T cell, B cell, NK, and myeloid-derived subpopulations. On average, we processed 2.12x10⁵ live cells per tissue sample (range = 3.73x10⁴ – 4.54x10⁵), with approximately 40% of events recorded identified as CD45+ (mean = 8.2x10⁴, range = 5.17x10³ – 2.43x10⁵). Using defined lists of immune cell lineage markers (Table 1), we first identified a total of 20 cell clusters composed of both myeloid and lymphoid lineages, in which we identified two CD3+CD8+ clusters, two CD3+CD4+ clusters, two CD3-CD19+ clusters, three CD16+CD56+CD3-CD19- clusters, three CD3+TCRgd+ T cell clusters, and six CD11b+ clusters (Supplementary Figure 2). Clusters with expression patterns associated with each primary lineage (CD4+ T cells, CD8+ T cells, CD11b+ myeloid cells, CD19+ B cells) were merged and re-clustered using markers specific to that cell type (Table 1; Supplementary Figures 3–6). Overall, this resulted in the identification of 66 distinct immune cell clusters across patients and tissues (Figure 1A).

We next sought to identify cell subsets that were enriched in patient tumors. To do this, we tested for paired differences in cell proportions between tissues within patients. We identified 57 cell subsets with different proportions between PBMC and tumor tissues (paired t-test with Bonferroni correction, q<0.05), 31 of which were specifically elevated in tumor relative to PBMC (Figure 1B). Of these, 14 were not different when comparing tumor to adjacent tissue indicating that their increased proportions relative to PBMC likely reflected recruitment to lung tissue generally and not to lung tumors specifically. We identified 39 subsets with significantly different proportions between matched tumor and adjacent tissue (paired t-test; Bonferroni, q<0.05), the majority of which (n=22, 53.8%) were elevated in tumor (Figure 1C). Notably, of those 21 subsets, 16 (76.2%) were among those found to also be elevated in tumor relative to PBMC. Additionally, 12 subsets were found to be significantly depleted in tumor on average, relative to both adjacent tissue and PBMC (Figures 1B, C). These cells included senescent CD8+ cells (CD45RO-CD45RA-CCR7-CD27-), as well as 2 NK subsets and 4 myeloid cells (classical monocytes, dysfunctional pDCs, cDCs, and granulocyte precursor cells).

We next focused on the 16 subsets whose proportions were elevated in tumor relative to both tumor-adjacent tissue and to PBMC, a pattern consistent with tumor infiltration. Notably, a majority, 10 (62.5%), of these tumor-infiltrating subsets were B lymphocytes corresponding to a variety of differentiated subtypes. These included resident activated memory B cells (CD19+CD20+IgD-CD21+CD27+CD69+), IgD+ activated memory B cells (CD19+CD20+IgD+CD21+CD27+CD69-), antigen presenting B cells (CD19+CD20+IgD-CD27-CD38+CD21+CD86+PD-1+), activated B cells (CD19+CD20+IgD-CD21+CD69-CD11c+CD86loCD79Blo), activated memory B cells (CD19+CD20+IgD-CD27+CD38-CD21lo), atypical resident B cells (CD19+CD20+IgD-CD27-CD38-CD21loCD69+), resident memory (CD19+CD20+IgD-CD27+CD38-CD21-CD69+), plasma cells (CD19+CD20-IgD-CD27+CD38++CD21-), IgD+ plasma cells (CD19+CD20-IgD+CD27+CD38++CD21-), and terminally differentiated plasma cells (CD19loCD20-IgD-CD27+CD21-CD38++). The remaining 6 subsets included CD4+CD69+ effector memory T cells (CD45RO+CCR7-CD27+CD69+), CD4+ CD69+ central memory (CD45RO+CCR7+CD27+CD69+), two regulatory CD4+ T cell subsets (CD45RO+CCR7-CD27+CD69+CD25+), CD8+ resident memory T cells (CD45RO+CCR7-CD27+CD69+) and CD8+ CD69+ TEMRA cells (CD45RA+CD45RO-CCR7-CD27+).

PCA using the proportions of the 66 subsets revealed distinct profiles for each tissue (Figure 1D). Based on principal components (PCs) 1 and 2 there was partial overlap between tumor and normal adjacent tissue clusters; however, PBMC samples formed a distinct cluster along PC1, which explained 34.7% of variation across samples. In contrast, adjacent and tumor tissues separated along PC2, which only explained 14.1% of variation across samples. PCA also revealed that the average inter-patient distances when comparing PBMC were much lower than those for either tumor or tumor-adjacent tissues (Figure 1E). Variability within tumor samples was especially pronounced for the 16 cell subsets whose enrichment patterns indicated they were tumor-infiltrating (Supplementary Figure 7). Hence, background interindividual variability seems likely to explain only a small part of the variability found in TMEs.

Our observation that tumor immune profiles are highly variable between patients (Figures 1D, E; Supplementary Figure 7) is consistent with previous studies that have shown variation in the numbers of tumor-infiltrating and tertiary lymphoid structure-associated immune cells among LUAD patients (8, 11, 19, 36, 37). Notably, many previous studies have noted the prognostic value of B cells in LUAD (12, 19, 38), and their potential interactions with other immune cell types, including CD4+ and CD8+ subsets (8, 11, 36). However, for the majority of these studies, the lower resolution profiling of immune sub-populations, particularly B cells subsets, has led to demarcation of a limited number of patient groupings with inconsistent prognostic value.

Here we sought to identify consistent tumor-infiltrating immune cell profiles among groups of patients, by performing unsupervised k-means clustering of our patients using tumor tissue proportions of the 66 CD45+ immune cell subsets described above. Using this approach, we identified 4 groups of patients (Groups 1-4), each characterized by distinct TI immune cell profiles. Using PCA, variation among patients and their respective assigned clusters was largely explained by PC1 (26.99%), PC2 (20.44%), and PC3 (15.34%; Figures 2A, B). We next used linear regression to determine which cell subsets explained the most variation in a given patient group relative to all other groups combined; this allowed us to effectively identify cell subsets that were either significantly enriched or depleted in each group. In total, we identified 35 immune cell subsets that wre significantly enriched or depleted (Bonferroni, q<0.05) in at least one group (Figure 2C; Supplementary Figure 8). In addition, to ensure these cell subsets were indicative of patient subgroup rather than other clinical features of these patients, we confirmed the lack of association with patient stage, biological sex, and age (Supplementary Figure 9; Table 3).

Patients in Group 1 had significantly higher proportions of 10 B cell subsets, including resident activated memory B cells (beta coeff. = 6.66; q=4.13^-5), naïve B cells (beta coeff. = 1.09; q=0.005), and antigen presenting B cells (beta coeff. = 0.272; q=0.013). Notably, 7 of these B cell subsets were among the 16 originally found to be enriched in tumors when we collectively analyzed all patients in this cohort (Figure 1B). The increased proportion of these B cell subsets in Group 1 was also reflected by the increased proportions of total B cells (CD19+ cells) relative to all other groups; indeed, in some Group 1 patients, CD19+ proportions represented >40% of all tumor-infiltrating CD45+ cells (Figure 2D). In addition, these increased B cell subsets were also associated with significantly higher proportions of CD4+ CD69+ central memory T cells (beta coeff. = 2.04; q=0.010) in Group 1 (Figure 2C). This combination of TI immune cell subsets, particularly the elevated representation of B cell populations was consistent with previous reports (8, 11). It is these particular B cell subsets that, when found enriched in tertiary lymphoid structures, have been associated with improved patient outcomes (8, 11, 19, 36). In contrast to Group 1 enriched cell subsets, CD4+ activated effector memory (beta coeff.=-0.270; q=0.024), classical monocytes (beta coeff.=-5.23; q=0.005), and mature NK cells (beta coeff.=-1.57; q=0.043) were moderately lower in Group 1 relative to the other groups combined; however, the proportions of these subsets in Group 1 were not uniquely lower when compared individually to Groups 2 and 3 (Supplementary Figure 10).

For Group 2 patients, tumors were significantly enriched for CD8+ resident memory T cells (beta coeff. = 13.0; q=4.22^-7), whereas the tumors of Group 3 patients were found to be uniquely enriched for CD4+ activated effector (Th1-like; beta coeff. = 25.6; q=9.95^-13) (Figure 2C). In both instances, these cell subsets were dominant subsets among tumors of Groups 2 and 3, respectively (Figure 2E). For example, in all three individuals of Group 3, CD4+ activated effector (Th1-like) cells represented >20% of all CD45+ subsets, compared to the other three patient groups in which the mean was <10% (Figure 2E). Finally, tumor profiles of Group 4 patients were dominated by myeloid cell populations, including classical monocytes (beta coeff. = 8.79; q=4.73^-9) and early myeloid-derived suppressor cells (early MDSCs) (beta coeff. = 1.36; q=0.018), as well as multiple NK cell populations, including mature NK cells (beta coeff. = 2.86; q=1.19^-5). In addition, many T lymphocyte populations were elevated in Group 4 relative to groups 1, 2, and 3; however, among these cell types there was biased representation of naïve and TEMRA subsets. These included naive CD4+ (beta coeff. = 0.294; q=0.035) and CD8+ (beta coeff. = 0.224; q=0.019) T cells, and senescent CD8+ T cells (beta coeff. = 1.61; q=0.003). Collectively, this signature was suggestive of an immunologically inactive tumor, consistent with earlier studies, potentially associated with the formation of myeloid-derived suppressor cells (43).

The most notable signature differentiating Group 1 from Groups 2–4 was the overall enrichment of CD19+ B cells among CD45+ cells (Figure 2D). This observation was consistent with previous reports showing elevated CD19/20+ B cells in some LUAD patients (8, 11, 19, 36). However, given the increased subset resolution possible with our dataset, we sought to further define the B cell populations underlying this enrichment. Specifically, we asked whether in addition to the enrichment of B cells in the tumor overall among Group 1 patients, whether the distribution of subsets within CD19+ cells was also distinct in Group 1. We reasoned that within-B cell subpopulation differences might explain why Group 2 and Group 3 show higher representation of CD8+ and CD4+ T cell subsets, further supporting differential immune cell interactions across groups.

Of the 15 B cell clusters identified, we identified 5 that showed statistically significant variation between the four patient groups. These included, resident activated memory B cells (ANOVA, p = 0.003), IgD+ memory B cells (ANOVA, p = 0.007), activated naive B cells (ANOVA, p = 0.009), activated B cells (ANOVA, p = 0.019), and antigen presenting B cells (ANOVA, p = 0.03) (Figure 3). Group 1 and Group 2 patients exhibit similar proportions across most B cell subsets, with a significant difference observed only in antigen presenting B cells (student’s t-test p<0.01; Figure 3). In contrast, Group 3 patients show notably lower proportions in four B cell subsets (antigen presenting, activated naive, atypical resident, and naive B cells) compared to Group 1 (student’s t-test p<0.05; Figure 3). Group 4 stands out as the most divergent, with two subsets (resident activated memory and antigen presenting B cells) present at lower proportions, and three subsets (IgD+ memory, activated naive, and activated) at significantly higher proportions relative to Group 1 (unpaired t-test p<0.05; Figure 3). Together these results demonstrate that not only are B cells elevated in particular patient groups, but that within the B cell compartment, there are specific distributions of B cell subpopulations.

To extend our investigation of B cell mediated immune responses in the LUAD TME, we next sought to characterize tumor-associated signatures in the expressed BCR repertoire. We conducted bulk AIRR-seq on IgM and IgG repertoires of matched PBMC, adjacent tissue, and tumor tissue from a subset of patients from the cohort studied above (n=29; Supplementary Table 3); this included 12, 7, and 4 individuals from Groups 1, 2, and 4 respectively. Given the infiltration of B cells noted in tumors (Figure 1), we reasoned that there would be associated changes in the BCR repertoire providing evidence of B cell activation and an antigen-driven immune response. We assessed this by comparing measures of clonal expansion and SHM between tissues. Using the repertoire polarization metric (39), defined as the number of clones that represent 80% of the total sequenced repertoire, we found that compared to PBMC, both adjacent tissue and tumor BCR repertoires were significantly more polarized (Figure 4A). This was true for both the combined IgM/IgG repertoire, as well as for IgM and IgG repertoires when considered separately (Figure 4A). Supporting this observation, we noted that repertoire diversity (Simpson’s diversity index) was also lower in adjacent tissue and tumor relative to PBMC (Supplementary Figure 11A).

We next assessed SHM frequency in the IgM and IgG repertoires across tissues. As expected, for all tissues we observed higher frequencies of SHM in the IgG repertoire relative to IgM (Supplementary Figure 11B). However, within each isotype-specific repertoire, we noted additional differences in SHM frequencies between tissues. For IgM, SHM frequencies in PBMC were low (mean = 0.004), likely reflecting the fact that the IgM repertoire captured in the periphery was largely representative of antigen-naive B cells. In contrast, SHM in both the adjacent and tumor samples was significantly higher (Figure 4B). Inter-tissue differences in IgG SHM frequencies were more subtle. However, SHM frequencies were significantly higher in adjacent tissue, relative to both PBMC and tumor repertoires (Figure 4B). These trends were consistent when we considered silent and replacement mutations separately; however, we did note a slight increase in the percent of replacement mutations (relative to silent mutations) in IgM and IgG for both PBMC and tumor repertoires compared to adjacent lung tissue (Supplementary Figures 11C, D). Collectively, these results provide evidence that tumor infiltrating B cells have undergone clonal expansion and potentially affinity maturation, suggestive of antigen-associated immune responses.

To further examine differences between BCR repertoires in the tumor relative to those circulating in the periphery, we evaluated biases in heavy chain BCR gene usage in tumor IgM and IgG repertoires compared to PBMC IgM. We reasoned that genes enriched in tumor repertoires relative to those in the periphery could represent BCRs engaging tumor antigen. We calculated usage frequencies of IGHV, IGHD, and IGHJ genes in PBMC IgM, as well as tumor IgM and IgG separately. To first broadly compare IgM gene usage between tissues, we conducted PCA of IGHV, IGHJ and IGHD gene frequencies, in matched IgM PBMC and tumor samples (Figure 5A). This revealed that tumor and PBMC IgM repertoires had broadly differentiated gene usage profiles (Figure 5A). Specifically, nearest neighbor distances between repertoires within tissues were much lower than matched distances within donors between tissues (Figure 5B). We next sought to identify specific genes enriched in tumor IgM repertoires relative to PBMC IgM representing gene usage shifts common across donors. Based on a paired t-test, we identified 2 IGHV genes (IGHV3-7; IGHV3-74) and 1 IGHJ gene (IGHJ4) significantly enriched in tumor IgM after Bonferroni correction (Figure 5C). Additional IGHV, IGHD, and IGHJ genes were significantly depleted (Supplementary Figure 12A).

We next compared PBMC IgM with tumor IgG repertoires, and PCA again demonstrated separation between the two tissues (Figure 6A). Nearest neighbor distances within tissues and pairwise distances within donors between tissues further supported that tumor IgG repertoires had diverged from PBMC IgM baselines (Figure 6B). Unlike tumor IgM, no IGHV genes were significantly enriched in tumor IgG relative to PBMC IgM, consistent with the notion that gene-usage shifts in tumor were not uniform across donors. However, five IGHD genes and one IGHJ gene exhibited significantly higher usage in tumor IgG, indicating some common patterns of change (Figure 6C). As with IgM, several IGHV, IGHD and IGHJ genes were significantly depleted in tumor IgG repertoires (Supplementary Figure 12B).

While the genes highlighted in Figures 5 and 6 represent those that exhibited average biases in tumor repertoires across donors collectively, this analysis did not account for genes that were commonly enriched in tumors of smaller subsets of patients. To assess this, we evaluated each patient individually, identifying the IGHV genes with the highest delta values between tumor IgM and IgG (separately) relative to PBMC IgM. In tumor IgM repertoires, 8 IGHV genes ranked among the top 5 most differentially used genes in at least 5 patients; the most frequently observed gene was IGHV3-7, identified in 19/26 patients (Figure 7A). Notably, this gene was also found to be significant overall in the matched comparison of PBMC and tumor IgM (Figure 5C). In tumor IgG, there was a greater number of IGHV genes (n=13) found among the highly ranked differentially used genes in 5 or more patients (Figure 7B), indicating again that the IgG tumor IGHV repertoire showed evidence of IGHV gene usage convergence only across smaller subsets of patients, rather than the majority of patients.

To test whether tumor-enriched IGHV genes in each patient coincided with changes in additional signatures associated with antigen selection, we assessed differences in: 1) the size of clones represented by those genes compared between tumor IgM/IgG and PBMC IgM; 2) SHM patterns; and 3) D and J gene usage profiles. Examples provided for the top two most commonly enriched IGHV genes in both tumor IgM (Figure 7C; Supplementary Figure 13A) and IgG (Figure 7D; Supplementary Figure 13B) demonstrated that clones associated with these genes are larger, have higher SHM frequency, and at an individual level, exhibit biased D and J usage in tumor relative to PBMC. This indicated that genes enriched in tumor reflected recruitment of particular clones composed of specific rearrangements of V, D, and J gene segments that had specifically expanded and taken on SHM. Notably, these trends were observed when we collectively considered features of BCRs associated with all top enriched genes relative to their depleted counterparts within the same donors (Figures 8A, B). Similar to gene-level data shown in Figures 7C, D, we observed that tumor-enriched IgM BCRs were represented by larger clones and increased SHM (Figure 8A). The same was also noted for IgG BCRs with respect to clone size differences (Figure 8B). However, we noted differential SHM signatures were specific to increased replacement (‘R’) mutations (Figure 8B), rather than overall SHM. We also tested for additional differences in CDR3 characteristics, including length, overall charge, and the fraction of basic or acidic positions. For IgM, we specifically noted that CDR3 length was shorter on average in enriched BCRs, but that CDR3 sequences in these BCRs also had a greater proportion of basic amino acids (Figure 8A). Likewise, enriched IgG BCRs were also associated with differential CDR3 features, relative to depleted BCRs; however, these features were different from those observed for IgM (Figure 8B).

Together, the results above demonstrate tissue-specific differences between tissues in the BCR repertoire, including broad trends that are noted across patients. However, given the distinct immune cell profiles observed in each of the identified patient groups (Figure 2), we reasoned that BCR signatures would also vary by group. To assess this, we first evaluated differences in BCR clonal structure between patient groups by comparing repertoire polarization of IgM and IgG in PBMC and tumor. In both IgM and IgG, Group 1 and Group 2 tumor repertoires were more polarized than Group 4 (ANOVA: IgM p=0.043; IgG p=0.005; Figures 9A, B). Notably, we observed that, across samples from each group, the top 10% of clones (by size) in tumor were predominantly represented by “enriched” genes, indicating that as expected these tumor-enriched genes were likely driving the group differences in the degree of clonal expansion (Figures 9C, D). Particularly within the IgG compartment, enriched genes made up over 99% of the largest clones in most group 1 (n=7/11) and 2 (n=5/7) patients (Figure 9D). However, despite differences in clonal expansion between groups, we did not observe significant differences in SHM profiles among this set of clones in each group (Figures 9E, F).

Our understanding of IgG subisotype variation in LUAD remains underexplored; however, it is well established that subisotypes can direct differential Ab-mediated immune responses. To expand on observations of bulk IgG repertoire profile differences between patient groups, we next leveraged a second higher-resolution AIRR-seq method that allowed us to characterize BCR repertoire features for each IgG subisotype in a subset of patients spanning Groups 1, 2, and 4 (n=18). This first revealed that Group 1 repertoires on average showed biased representation of IGHG1 and IGHG4 transcripts (Figures 10A–C), compared to Groups 2 and 4, which presented modest biases of IGHG2 and IGHG3. Additionally, we observed that a subset of patients within Group 1 uniquely showed elevated proportions of IGHG4 transcripts, in one instance exceeding 40% of the IgG repertoire (Figure 10B). Additionally, while we noted that tumor-enriched genes were found to represent large fractions of each subisotype, there were notable between-group differences in the representation of these genes among IGHG4 transcripts, revealing a particular elevation in Group 1 patients (Figure 10D). This was regardless of whether overall IGHG4 usage was elevated in a given patient. Notably, in two Group 1 patients, tumor-enriched IGHV genes comprised >75% of IGHG4 transcripts, dominated by a single gene that was uniquely elevated only in the IGHG4 compartment (Figures 10E, F).

These results suggested that a subset of B cells in these patients had undergone class-switching to IGHG4 from IGHG1. Our use of FLAIRR-seq data allows us to gain resolution of IGHC gene allelic variation within the IgG repertoire, and thus link expressed VDJs to specific IGHC alleles. In the case of Patient 915, who was heterozygous for both IGHG1 and IGHG4 (Supplementary Figure 14A), we were able to link the dominant IGHV gene (IGHV1-24) to single IGHG1 and IGHG4 alleles determined to reside on the same chromosome (“Haplotype 1”; Figure 11A; Supplementary Figure 14B), providing direct evidence of class-switch recombination. We observed that the proportion of IGHG4-IGHV1–24 transcripts derived from this haplotype were 4.6-fold more abundant than those using the IGHG1 allele from the same haplotype. This was in contrast to IGHV1–24 transcripts expressed from Haplotype 2, which were far less abundant, and equally distributed between IGHG1 and IGHG4 (Figure 11A). IGHV1–24 transcripts from Haplotype 1 were represented by 51 clones, but dominated by three clones (Figure 11B), two of which shared the same IGHJ gene (IGHJ4), but distinct IGHD genes. In these three clones, IGHG4 transcripts were observed at a greater proportion, and in all cases showed higher amino acid replacement mutations than IGHG1 transcripts, and higher rates in CDRs than in FWRs (Figures 11C, D). Together, this analysis provides evidence that tumor-enriched genes in Patient 915 are overrepresented among IGHG4 transcripts, resulting from class-switch recombination and IGHG4-specific clonal expansion, and associated with increased SHM.

Although increased tumor-derived IgG4 is generally linked to immunosuppressive environments and poor prognosis in other solid cancers (40–42), we next asked whether IGHG4 frequency in Group 1 patients correlated with specific TI immune cell profiles. However, one study in LUAD found IgG4 could be positively associated with outcome depending on driver mutations, possibly due to its inability to form immune complexes (21). Other studies in LUAD have either not distinguished IgG subisotypes or omitted IgG4 altogether because of low abundance (10, 12). Here, we next asked whether IGHG4 frequency in Group 1 patients correlated with specific TI immune cell profiles. In our cohort, IGHG4 frequencies were positively correlated with four TI subsets: cytotoxic NK cells (R = 0.83, p = 0.059), antigen-presenting B cells (R = 0.78, p = 0.013), senescent CD8+ T cells (R = 0.74, p = 0.022), and effector memory TCRgd cells (R = 0.68, p = 0.043; Figures 12A, B). These results suggest that elevated IGHG4 may reflect coordinated changes in cytotoxic and antigen-presenting populations driven by local class-switching and clonal expansion, linking IGHG4 expression to shaping of the TI immune microenvironment. Notably, although not statistically significant, Group 1 patients with high IGHG4 shared a trend toward poorer overall survival when compared to other patient groups used in FLAIRR-seq as well as all other patient groups (Supplementary Figures 15A, B).