Comprehensive multiomic characterization of human papillomavirus-driven recurrent respiratory papillomatosis reveals distinct molecular subtypes

Abstract

Recurrent respiratory papillomatosis (RRP) is a debilitating neoplastic disorder of the upper aerodigestive tract caused by chronic infection with low-risk human papillomavirus types 6 or 11. Patients with severe RRP can require hundreds of lifetime surgeries to control their disease and pulmonary papillomatosis can be fatal. Here we report the comprehensive genomic and transcriptomic characterization of respiratory papillomas. We discovered and characterized distinct subtypes with transcriptional resemblance to either a basal or differentiated cell state that associate with disease aggressiveness and differ in key molecular, immune and APOBEC mutagenesis profiles. Through integrated comparison with high-risk HPV-associated head and neck squamous cell carcinoma, our analysis revealed divergent molecular and immune papilloma subtypes that form independent of underlying genomic alterations. Cumulatively our results support the development of dysregulated cellular proliferation and suppressed anti-viral immunity through distinct programs of squamous cell differentiation and associated expression of low-risk HPV genes. These analyses provide insight into the pathogenesis of respiratory papillomas and provide a foundation for the development of therapeutic strategies.

Fig. 1. Copy number alterations are rare events in RRP.

Line graphs show the fraction of samples with indicated copy number alterations at corresponding genomic locations in a HPV-associated HNSCC, b RRP and c normal mucosa. Red and blue lines correspond to amplifications and deletions, respectively. CN = 0 corresponds to a complete loss of both chromosomes. HPV, human papillomavirus; HNSCC, head and neck squamous cell carcinoma; RRP, recurrent respiratory papillomatosis.

Fig. 2. HNSCC-associated mutations are less frequent in RRP.

Heatmaps shows mutations within 30 COSMIC CGC genes (rows) most frequently mutated in HPV-associated HNSCC in a normal mucosa (n = 7) b RRP (n = 21) and c HPV-associated HNSCC (n = 74). The barplots above each heatmap show the total number of mutations within each sample. Source data used to generate this figure are presented in Supplementary Table 1 and Supplementary Data 1. COSMIC CGC, catalog of somatic mutations in cancer—cancer gene census.

Fig. 3. Analysis of genetic and clinical characteristics reveal HPV type as predictor for development of pulmonary disease.

Heatmap shows mutations in COSMIC CGC genes (rows) in RRP samples (n = 21) (columns) predicted to be possibly or probably damaging based on PolyPhen score (REF). The bar plots on top show the number of clinically indicated interventions in the twelve months prior to biopsy (surgery frequency) and the total number of lifetime surgeries for each patient. Pulmonary disease status, HPV type and age of disease onset for each patient are shown below the heatmap.

Fig. 4. HPV type-specific gene expression in RRP and HNSCC.

Heatmap shows individual HPV gene expression of HPV types 6, 11, 16, 33 and 35 (rows) in RRP (n = 15), HPV-associated HNSCC (n = 55) and normal mucosa (n = 4) (columns). Color represent log2(TPM + 1). The top scatter plot shows the number of HPV integration events for each sample as determined by HPV gene transcript fusions. Only HPV-associated HNSCC samples with added HPV gene expression of 10 TPM or greater are shown. TPM, transcripts per million.

Fig. 5. Transcriptional variability in RRP associated with HPV type and intervention frequency.

a Scatter plot shows projections of samples onto the first and second principal components obtained by PCA. Shapes represent HPV type and color corresponds to the number of clinically indicated interventions in the twelve months prior to biopsy (surgery frequency). b Volcano plot shows results of differential gene expression analysis between PC1-low (PC1 projections < 0) and PC1-high samples (PC1 projections > 0) in a. Genes with an absolute log2 fold expression change greater than or equal to one and an adjusted p-value less than or equal to 10^-5 were considered to be differentially expressed as indicated by vertical and horizontal dashed lines, respectively. Using these criteria 561 and 1128 genes were up- and down-regulated in PC1-high samples, respectively. Genes with mean TPM ≥ 8 were considered for the PCA. c Dotplot shows Reactome pathway terms enriched in genes differentially expressed in either PC1-low (low frequency subtype) or -high (high-frequency subtype) samples. d Heatmap shows row-standardized expression of genes related to stemness identity, cell cycle progression and epithelial differentiation (rows) within low- (basal subtype) and high- frequency (differentiated subtype) of intervention samples as indicated in the top annotation. HPV type for each sample is shown above the heatmap. PCA, principal component analysis.

Fig. 6. HPV gene expression is associated with increased intervention frequency.

a Scatter plot shows expression changes of HPV genes comparing basal- and differentiated subtype samples. Log2 fold changes greater than 0 (dashed vertical line) indicate increased expression in differentiated subtype. P-values signifying differential expression are shown on the vertical axis. HPVx corresponds to HPV6 or HPV11. b Boxplot shows the expression of HPV genes within RRP samples stratified based on the number of clinically indicated interventions in the twelve months prior to biopsy, as indicated by color. HPVx corresponds to HPV6 or HPV11. P-value = 3.3 × 10⁻¹⁴, two-way ANOVA. Source data for this figure are presented in Supplementary Table 1 and Supplementary Data 2.

Fig. 7. Divergent immune function within RRP molecular subtypes.

a Heatmap shows row-standardized expression of genes related to interferon signaling, antigen processing, antigen presentation, T cell function as well as APOBEC and HPV gene expression (rows) within basal and differentiated subtype samples. HPV type for each sample is shown above the heatmap. b Boxplots show the fraction of CD4-positive cells (left panel) and CD8-positive cells (right panel) with the basal and differentiated subtype samples (p values are based on Wilcoxon rank-sum test). c The representative photomicrographs show IHC staining of CD4-positive cells (first column) and CD8-positive cells (second column) in a basal (first row) and differentiated subtype sample (second row). Each image shown is photographed at 200X magnification. Source data for this figure are presented in Supplementary Table 1 and Supplementary Data 4. IFN, interferon; APOBEC, apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like; IHC, immunohistochemistry.