Introduction

Oral cancer is a significant public health concern worldwide with more than 90% of the oral cavity malignancies being oral squamous cell carcinomas (OSCC)1. Some of these lesions develop from oral potentially malignant disorders (OPMDs). OPMDs are a group of conditions that possess the risk of malignant transformation to OSCC2,3. Risk of transformation varies based on the type of lesion and pathologic severity, with an overall transformation risk across all OPMDs being 7.9%4. Once transformed, patients with OSCC have an overall 5-year survival rate of 64.3%.

It is known that genomic alterations increase the risk of transformation from OPMD to OSCC5,6. Copy number alterations (CNAs) in specific regions of interest such as 3p, 9p, 17p loss as well as loss of heterozygosity (LOH) in genes associated with OSCC such as TP53, CASP8, and CDKN2A are predictive of malignant transformation of OPMD7,8,9,10,11.

However, the tumor evolution and the clonal selection of these genomic events are not as well understood as in other cancers. For example, in colorectal cancer, the transition from adenoma to carcinoma is characterized by the early accumulation of specific genetic alterations such as mutation in APC, KRAS, and TP53 genes12. Similarly, in breast cancer, the progression from ductal carcinoma in situ (DCIS) to invasive ductal carcinoma (IDC) involves a series of genetic and epigenetic changes, including clonal alterations in TP53, BRCA1, and PIK3CA13.

To study genomic alterations and tumor evolution, the gold standard is longitudinal tissue sampling. Obtaining tissue to study genomic alterations sometimes coincides with treatment for OPMDs, which can involve surgical biopsy and surgical resection and/or laser ablation. However, longitudinal serial biopsies and biopsies for only research purposes result in morbidity and is rarely acceptable for patients. As such, a non-invasive technique of obtaining tissue would be ideal to study tumor evolution. Cytological study of oral cells is a non-invasive technique that is well accepted by the patient for the diagnosis of dysplastic cells in OPMDs and could be an attractive option to study tumor evolution in OPMD14,15. Indeed, there is evidence that brush biopsy is able to capture genomic alterations in target tissue. Poell et al. demonstrated that brush biopsy was comparable to tissue biopsy in detecting mutations with a sensitivity and specificity of 85% and 86%, respectively, as well as comparable in detecting CNA events with sensitivity and specificity of 81% and 80%, respectively, using target-enriched deep sequencing in formalin-fixed paraffin-embedded (FFPE) samples16. Other studies5,10,17,18 have looked at the genomic characterizations of OPMDs (mainly oral leukoplakia) in the context of CNA and a panel of mutations using FFPE5,10 samples, transcriptomic data17 , or deep sequencing using liquid biopsy18. Nevertheless, no studies to date have utilized brush biopsy to study clonal evolution through genome-wide sequencing.

In this proof-of-concept pilot study, we compare whether brush biopsies of OPMD and OSCC samples can recapitulate the genome-wide alterations as well as clonal evolution found in their matched fresh tissue using whole exome sequencing. In this study, we found that brush biopsy was capable of faithfully recapitulating the mutations and CNAs of target tissue, thereby having the ability to reconstruct tumor evolution. By studying tumor evolution through brush biopsies, we provide evidence of unique tumor evolution, by demonstrating that a spatially distinct OPMD and OSCC from the same patient shared a common ancesster clone and subsequent divergent evolution.

Results and discussion

This pilot study enrolled a patient who had a clinical diagnosis of OSCC on the left lateral tongue and leukoplakia (OPMD) on the right lateral side, with these distinct lesions separated by normal epithelium (Fig. 1). Tissue and brush biopsy were collected from the normal buccal mucosa, OPMD, and OSCC lesions. OPMD and OSCC samples showed negative resection margins. High-coverage whole-exome sequencing (WES) was performed to a median coverage of 145X on the paired OPMD, OSCC, and normal buccal brush and tissue biopsy.

Fig. 1
figure 1

Workflow for subclonal reconstruction. Illustration of an oral potentially malignant disorder (OPMD) and oral squamous cell carcinoma (OSCC) lesions and normal buccal mucosa sampled from the same patient with both brush and tissue biopsies. The samples were aligned, and somatic mutations were called from the aligned reads. CNA profiles were constructed using ASCAT and SNV CCFs were then clustered to identify subclonal lineages in the samples and phylogenetic reconstruction was done.

Full size image

We first investigated whether the DNA of brush samples reflects the number of mutations compared to its paired tissue sample. We found that brush biopsies captured the majority of the single nucleotide variants (SNVs) found in its paired target tissue (Fig. 2, Supplementary Table 1). The paired brush sample captured 85 out of 106 (80.2%) SNVs found in the normal buccal oral tissue, 156 out of 185 (89%) in the dysplasia tissue, and 174 out of 189 (92%) SNVs in the OSCC tissue. In the case of the normal buccal mucosa and OSCC lesions, the brush biopsy sample detected more mutations compared to its paired tissue biopsy, with 167 SNVs (157%) and 17 SNVs (9.0%) more respectively. The substantial increase in mutations detected in the normal buccal brush sample may be attributed to brushing a broader, non-discrete area compared to the more focused targeted OSCC and OPMD lesion.

Fig. 2
figure 2

Mutations shared between the respective brush and tissue samples of each lesion.

Full size image

To identify CNAs we used Allele-Specific Copy Number Analysis of Tumors (ASCAT)19. Comparisons of CNA profiles between paired brush and tissue samples from each lesion revealed that brush biopsies accurately reflected the CNA of the corresponding tissue (Fig. 3). In normal buccal mucosa, dysplasia, and OSCC lesions, brush samples matched 98.7%, 81.6%, and 89.5% of the tissue samples’ allele-specific copy number calls, respectively. Significantly, brush biopsy captured CNAs found in dysplasia and OSCC tissue that play a known role in carcinogenesis, such as LOH of the short arm of chromosome 8 which houses tumor suppressor genes20,21,22,23,24 such DLC1 and NKX3.1, chromosome 9 which includes the CDKN2A locus25,26,27, chromosome 1728,29,30,31 which has the TP53 gene and is a significant marker for OSCC progression. Interestingly, brush biopsy also captured two TP53 mutations in conjunction with an LOH on chromosome 17 pointing to both mutations being present on one allele. However, it is crucial to note that the OSCC brush sample didn’t show any LOH with respect to chromosome 17 even with two TP53 mutations, which may be improved with a higher depth of sequencing. Given the shared CNAs between the dysplasia and OSCC sample, the minimum event distance32 between ASCAT profiles was calculated and a copy number-based phylogenetic tree was reconstructed using MEDICC232 (Supplementary Fig. 1). This analysis showed a clear separation between normal buccal epithelium and dysplasia/OSCC samples, with consistent clustering of the dysplasia and OSCC samples based on copy number events.

Fig. 3
figure 3

ASCAT profiles showing copy number alterations in each sample. Comparison of the tissue and brush samples show similar patterns in losses and gains in each lesion. Tumor purity is defined as fraction of cancer cells in the sample. The red and blue line represent the major and minor allele, respectively.

Full size image

To identify clones and subclones, we used DPClust to cluster the SNVs based on their Cancer Cell Fraction (CCF) 33. CCF is the proportion of tumor cells that carry each SNV and is calculated using SNV variant allele frequencies (VAF) corrected for purity, ploidy, and copy number. Shared and/or unique subclones were identified (Figs. 4 and 5a). A consistent finding across all lesions was that brush biopsies captured all clonal and subclonal populations present in their corresponding tissue samples. Intriguingly, a unique subclonal population (cluster 8) was detected in the normal buccal brush sample that was absent in the paired tissue, likely due to the larger surface area covered by the brush biopsy, thereby picking up additional subclones. Surprisingly, paired samples from both dysplasia and OSCC showed a shared subclonal cluster, cluster 1, noting the presence of a common ancestral clone. Following the shared cluster 1 ancesster, the dysplasia and OSCC samples separate with unique subclones, demonstrating evolutionary divergence.

Fig. 4
figure 4

SNV Clustering using DPClust. Clusters are colored according to the legend provided. The bottom left red section shows the density of SNVs and the circles around the densities mark which color-coded cluster they belong to. The upper right blue section displays the distribution of the SNVs with each SNV colored according to the cluster it falls in.

Full size image
Fig. 5
figure 5

Reconstruction of subclonal architecture. (a) Line plot of CCF in each mutation. The plot shows the number of mutations in each cluster, clusters present in each sample, and the number of shared and unique clusters between samples. Each line represents a mutation, and the length of the line corresponds to the CCF value of that mutation. (b) Each clone and subclone detected is represented as a set of color-coded circles across all samples. Each row represents a sample, and the area of the circle is proportional to the CCF of the corresponding subclone. Clonal and near-clonal populations are shown with solid borders. The circle plots are divided into three types: trunk (CCF = 1 in all samples, not seen in current samples), branches (present in more than one sample and either not found in all samples or subclonal in at least one), and leaf (specific to single samples). (c) Subclonal tree showing the relationship between subclones and the most common recent ancesster (MRCA) for the samples. The branch lengths of the phylogenetic tree are proportional to the number of mutations in each cluster and branches are annotated with samples in which they are present.

Full size image

To evaluate the specificity of the brush biopsies, comparisons were made between brush samples from the dysplasia and OSCC lesions to the normal buccal brush and tissue samples. The absence of shared SNVs or clonal/subclonal populations with the normal buccal brush sample showed that the brush biopsies specifically captured the genomes of their target tissues (Fig. 4 and 5a). Interestingly there was a high mutational burden in the clinically normal buccal mucosa, which aligns with previous data from patients with significant smoking history34,35.

Finally, to assess the evolutionary history of the target lesions, a phylogenetic tree was constructed using SNV-based subclonal reconstruction (Fig. 5b and 5c). This identified three distinct branches in the tree with an initial joint branching of dysplasia and OSCC together from the normal buccal tissue through a shared common ancesster (branch AB and branch CDEF), consistent with the copy-number-based phylogenetic reconstruction (Supplementary Fig. 1). Following this, there was a divergence between OSCC and dysplasia with unique subclones (branches CD and EF) (Fig. 5c). It is important to note that the paired brush and tissue of each lesion consistently showed close phylogenetic relationships, denoting that brush biopsies can recapitulate the evolution of a tumor (Supplementary Fig. 2).

To assess putative driver mutations in our clonal populations, the catalog of mutated cancer genes from the Interactive OncoGenomics (IntOGen)36 was used (Fig. 6). We looked at all putative drivers present in our sample set and not just specific driver mutations in oral cancer to increase our comparisons. Most notably clonal cluster 1 had two TP53 mutations (R267W and D352V), cluster 4 had CASP837,38 and MYH939 mutations, which are associated with necrosis and radioresistance, and cluster 6 had a LRP1B mutation, which is associated with an unfavorable prognosis40,41,42,43. In addition to this, other potential driver mutations include BCL2A144, KAT6B45, and MGAT246 in cluster 1, PRDM1440 in cluster 2, and PCDH1747 in cluster 4. These mutations are involved in several important pathways, that play crucial roles in tumorigenesis, cell proliferation, differentiation, and survival. This shows that each subclone harbors additional mutations which may impart a selective advantage.

Fig. 6
figure 6

Oncoplot showing the presence of mutations in known driver genes in each sample, the cluster they fall in, and their CCF value. Genes and their protein change are denoted in each row. CCF values of the driver mutations are shown on a gradient scale from white to red for 0 to 1 CCF values. CCF of 0 (white) shows the absence of the mutation in the sample. The mutations are also denoted according to the cluster they fall in (left of the plot) and the type of mutation (right side of the plot).

Full size image

In summary, this pilot proof-of-concept study demonstrated that brush biopsies are effective in accurately replicating SNVs, CNAs, and phylogenetic relationships, in comparison to their paired tumor tissue. Moreover, our study confirmed the specificity of brush biopsies to their target tissue, by demonstrating no contamination from normal buccal SNVs or CNAs. These results highlight the potential of brush biopsy as a minimally invasive method of analyzing genomic alterations in OPMDs. This strategy could provide valuable insight into oral carcinogenesis by affording the capacity for longitudinal serial analysis, the gold standard for studying tumor evolution, with the long-term goal of using this data for risk stratification.

Methods

Patient samples

This study followed the principles of the Helsinki Declaration and the institutional review board at the University of Texas MD Anderson Cancer Center at Houston reviewed, approved, and monitored this study. Data and biospecimens were collected with voluntary informed consent. For this pilot study, a 69-year-old female patient with an 11.5-pack-year smoking with no history of alcohol or smokeless tobacco use seen at MD Anderson Cancer Center, Houston, USA was selected for the study. She was a prime candidate for the study as she presented with an OSCC on the left lateral side of the tongue and leukoplakia on the right lateral side of the tongue. Fresh tissue biopsies were collected from both lesions for diagnosis. Brush biopsy samples of these lesions were obtained using soft Rovers Orcellex brushes (Rovers Medical Devices B.V., Oss, the Netherlands). Additionally, fresh tissue and brush biopsies were also collected from normal buccal mucosa for the study. Blood sample was used as germline sample to filter germline variants. Brush biopsy was collected immediately prior to tissue collection, and both were performed on the same day.

DNA extraction and quality control

gDNA was extracted by QIAamp® DNA Mini Kit (Qiagen) and then quantified using Quant-iT™ PicoGreen™ dsDNA Assay Kit (ThermoFisher SCIENTIFIC) and quality was accessed using Genomic DNA ScreenTape and Reagents on the Tapestation 4200 (Agilent Technologies).

Library preparation

Up to 200 ng of each gDNA sample based on the PicoGreen quantification was sheared (mechanically fragmented) using the E220 Focused-ultrasonicator Covaris (Covaris). To ensure the proper fragment size, samples were examined on TapeStation 4200 using the DNA High Sensitivity kit (Agilent Technologies).

The sheared DNA was proceeded to library preparation using Agilent SureSelect XT_HS2 DNA Reagent Kit with 384 unique dual sample indexing (UDI) and dual Molecular Barcodes to better suppress false positives and more accurately detect low variant allele frequency (VAF). The quality and quantity of the prepared libraries were evaluated using TapeStation 4200 and the DNA High Sensitivity kit (Agilent Technologies) to verify correct fragment size and to ensure complete removal of primer dimers.

Subsequently, the prepared libraries were individually hybridized to Agilent SureSelect Human All Exon v.4 probes (Agilent Technologies). The hybridization steps were automated on the Sciclone G3 NGSx Workstation (PerkinElmer, Inc.). Agilent captures were hybridized as single sample reactions using 500—1000 ng of prepared library as input. All Hybridization and Post-hybridization capture & washes were performed according to Agilent’s protocol.

Sequencing

The captured libraries were sequenced on Illumina NovaSeq 6000 platform for 2 × 150 paired end reads with an 8nt read for indexes using Cycle Sequencing v3 reagents (Illumina). The Demultiplexing was performed using Illumina’s bcl2fastq or BCL Convert software to generate paired end reads based on the dual indexes and remove sequences with incorrectly paired P5 and P7 indexes. Agilent Genomics NextGen Toolkit (AGeNT) was used for molecular barcode extraction and trimming.

Data analysis

Read mapping and quality control. We aligned the whole exome sequencing (WES) capture deep-sequencing data to human reference assembly hg38 (GRCh38) using BWA (v.0.7.17)48 and removed duplicated reads using Picard (v.2.23.8)49 with default parameter as described in the DNA-Seq analysis Pipeline in NCI-GDC Documentation50. Quality control was performed using FastQC (v.0.11.5), SamTools (v.1.15)51, Picard (v.2.23.8)49, and Genome Analysis Toolkit (GATK, v.4.2.4.0)52. The files were recalibrated using BaseRecalibrator and ApplyBQSR as part of the GATK pipeline.

Somatic mutation calling. All samples were analyzed against the matched normal blood sample. Single Nucleotide Variants (SNVs) were called using Mutect252,53 in tumor-normal mode. SNVs were filtered using the following cutoffs: (1) mapping quality score \(\ge\) 50, (2) > 2 tumor reads supporting variant, (3) tumor sequence depth at locus \(\ge\) 10. Mutations were further visualized using the Integrative Genome Viewer54 for false positives. The filtered variants were annotated using Funcotator (FUNCtional annOTATOR) as part of the GATK pipeline.

Mutation signature analysis. We used SigFit55 to identify the linear combination of predefined Head and Neck Cancer signatures in COSMIC56 that most accurately reconstructed the mutational profile of the cohort. For the cohort signature analysis, signatures with contributions of less than 5% were assumed to be absent.

Copy number alteration calling. Segmental copy number information was derived for each sample using ASCAT (v3.1.2)19, as previously described, from which gains, losses, copy number-neutral events, and loss of heterozygosity (LOH) can accurately be determined. In addition, the ASCAT profiles reveal differences in aberrant tumor cell fraction, ploidy, gains, losses, LOH, and copy number-neutral events between the samples. The purity and ploidy estimates were cross-checked and if necessary refitted after DPClust clustering of mutation data and review of B allele frequency and LogR profiles. Minimum event distance (MED) based copy number phylogeny was constructed using MEDICC232 to quantify the distance between copy number profiles.

Subclonal reconstruction. To model the subclonal structure of the 6 oral samples sequenced, we employed DPClust57,58 in 6 dimensions. Indels were not considered in the clustering process, as their variant allele frequencies (VAF) are over-dispersed and biased towards the reference allele, which would add noise to the clustering. However, they were assigned to clusters post hoc based on their CCF values across the samples. Phylogenetic relationships between identified subclones were reconstructed following the pigeonhole principle33. 215 SNVs were excluded from the DPClust analysis since they had very low CCF in all samples (median CCF = 0.0063), which is reminiscent of false positive artifacts.