Supplementary Figure S2. Antigen presentation module genes were enriched among genes negatively correlated with hsa-miR-149-5p in lung-related datasets.

Clinical annotations for endobronchial biopsies that were profiled by miRNA-seq stratified by PML molecular subtypes.

Supplementary Figure S10. Analysis of IMC data from the epithelial tissue.

hsa-miR-149-5p is highly expressed within the epithelium and associates with PML progression. A, Bubble plot showing the correlation between the residual expression level of hsa-miR-149-5p and cell type marker genes, including CD3G (T cells), CD19 (B cells), CD68 (macrophages), KRT5 (basal cells), FOXJ1 (ciliated cells), MUC5AC (goblet cells), SCGB1A1 (club cells), MUC5B (secretory cells), and CEACAM5 (peri-goblet cells). Correlation results for hsa-miR-34b-5p, hsa-miR-449c-5p, and hsa-miR-150-5p are shown as controls. B, Enrichment of hsa-miR-149-5p across cell type compartments in the FANTOM5 project. Bar chart shows the normalized expression levels of hsa-miR-149-5p in the FANTOM5 samples (n = 399) (top). Vertical bars indicate the position of FANTOM5 samples derived from the epithelial (red) or immune (orange) cell compartments (bottom). C, Scatter plot of the Pearson correlation between the normalized expression levels of hsa-miR-149-5p and NLRC5 within the FANTOM5 samples derived from the epithelial cell compartment (n = 87). The dashed red line represents the linear regression fit and the shaded region indicates the 95% confidence interval. D–F, The hsa-miR-149-5p density was detected by miRNA-ISH from proliferative-subtype PMLs (n = 19; 10 progressive/persistent and nine regressive). D, Boxplot showing the hsa-miR-149-5p density per area between the epithelium and nonepithelium regions. E, hsa-miR-149-5p density per nucleus between the progressive/persistent and regressive PMLs within regions of normal epithelium, regions with hyperplasia or metaplasia (hyper-metaplasia), and regions of bronchial dysplasia. F, Representative miRNA-ISH staining of hsa-miR-149-5p for progressive/persistent PMLs (top) and regressive PMLs (bottom). Arrows indicate stained hsa-miR-149-5p within epithelium. Boxplots indicate median with IQR and whiskers indicate minimum and maximum measurement in (D) and (E). P values were FDR (Benjamini and Hochberg) adjusted and determined by Pearson correlation coefficients (A), the two-sided paired Wilcoxon test (D), and a linear mixed-effects model (E). *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ****, P ≤ 0.0001; ns, no significance.

Study overview. A, High-quality miRNA-seq data was generated from 156 frozen endobronchial biopsies (n = 28 patients). These samples also had paired mRNA sequencing data. Among the 156 biopsies, 19 adjacent FFPE biopsies, transcriptionally defined as the proliferative molecular subtype, were sectioned and stained to assess spatial localization of key findings. B, The paired miRNA–mRNA sequencing data were used to construct a miRNA–gene network based on database evidence of miRNA target genes and significant negative correlation between miRNA and gene expression patterns (gray circles are miRNAs and red, green, and yellow circles are genes for which the color represents their membership in previously described gene coexpression modules that define PML molecular subtype). Target information was obtained by combining a sequence-based prediction database and experimentally validated databases. Pearson correlation was calculated from expression residuals for each miRNA–mRNA pair, and those with significant negative Pearson correlation coefficients (FDR ≤ 0.05) were selected. Next, we filtered edges to identify miRNAs that may have central gene expression regulatory roles within the miRNA–gene network in each gene module. For each gene module, we identified and retained only the connections from miRNAs for which (i) there are more significantly negative correlated target genes in the module than in other modules (Fisher exact test, odds ratio>1) and (ii) the Fisher-Z–transformed Pearson correlation coefficient densities with its targets within the module were more negative than its targets in other modules (t test, t < 0). miRNAs with significant results from both tests (FDR ≤ 0.05) were considered module-associated miRNAs, and other miRNAs were removed from the network. C–E, For each FFPE biopsy, three consecutive sections were stained for different purposes: (C) H&E to annotate the histologic grades present in the epithelial tissue, (D) chromogenic miRNA-ISH to quantify hsa-miR-149-5p expression and localization, and (E) IMC using 28 metal-tagged antibodies to identify epithelial and immune cell types and hsa-miR-149-5p target protein expression and localization. Single or multiple ROIs per sample were selected for IMC ablation in 14 samples. The epithelial histology in both miRNA-ISH and IMC images was annotated based on the corresponding adjacent H&E section. The schematic was constructed using data images from Fig. 4. [Created in BioRender. Chiu, D. (2026) https://BioRender.com/g78j632]

Supplementary Figure S11. Analysis of IMC data from the non-epithelial tissue.

Hsa-miR-149-5p and NLRC5 target genes are associated with PML progression, and in vitro modulation of NLRC5 changes MHC-1 expression and CD8 T cell–mediated cytotoxicity. A and B, Boxplots of residual expression levels of progressive/persistent (dark green) and regressive (light green) proliferative subtype PMLs from the discovery cohort (GSE109473) for (A) hsa-miR-149-5p significant negatively correlated target genes in module 9 and (B) NLRC5 target genes. C and D, SW900 cells were transfected with mock control (red orange), hsa-miR-149-5p (green), or siNLRC5 (blue) for 48 hours. Relative gene expression detected by qRT-PCR (6–10 replicates per condition) of (C) NLRC5 and its downstream targets (HLA-A, HLA-B, HLA-C, TAP1, and PSMB9) and (D) hsa-miR-149-5p. PPIA was used as the reference gene to calculate relative fold change values. E and F, Barplot showing the difference in (E) mean fluorescent intensity (MFI) of surface MHC-I and (F) percent cytotoxicity between IFNγ-simulated and unstimulated conditions between mock- and siNLRC5-treated SW900 cells. Data indicate median with IQR, and whiskers indicate minimum and maximum measurement. Error bars indicate standard errors. P values were determined using mixed-effect linear models (A, B, and E), two-sided Wilcoxon tests (C and D), and two-sided t test (F). *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001; ns, no significance.

Supplementary Figure S6. Correlation between hsa-miR-149-5p and NLRC5 expression levels within the samples from the FANTOM5 project.

Protein expression of NLRC5 and its targets is correlated, and basal cells with increased NLRC5 expression are in close spatial proximity to CD8 T cells. A, tSNE visualization of eight epithelial and four immune cell populations within the epithelium (n = 41,147 cells) from IMC data (n = 14 samples; seven progressive and seven regressive). B, Heatmap of the average log-Z–normalized expression within the epithelium of canonical markers (left) and the average log-Z stromal-normalized expression within the epithelium of NLRC5 and its downstream targets (right) in identified cell types from (A). Bar plot shows the relative proportion of each cell type. C, Bubble plot shows the differential composition of cell types in hyper-metaplasia and dysplasia as compared with normal epithelium. Color represents the compositional estimate, and dot size represents the log FDR value computed by sccomp. D, tSNE visualization of epithelial cell populations colored by NLRC5-high (dark blue) vs. NLRC5-low (light blue) expression. E, Violin plot of the log-Z stromal-normalized expression of NLRC5 target genes (TAP1, PSMB9, and B2M) stratified by the NLRC5-high and NLRC5-low labels. F, tSNE visualization of eight immune and three stromal cell population in the stroma region (n = 46,254 cells) from IMC data (n = 14 samples; seven progressive and seven regressive). G, Heatmap of the average log-Z–normalized expression within the stroma of canonical markers (left) and the average log-Z stromal-normalized expression within the epithelium of NLRC5 and its downstream targets (right) in identified cell types from (F). Bar plot shows the relative proportion of each cell type. H, Violin plot of log-Z stromal-normalized expression of NLRC5 in the epithelium (red) vs. the stroma or nonepithelial tissue (yellow). I and J, Cellular neighbors for each cell were defined based on a centroid distance within 30 μm. Cell–cell proximity score within each ROI (n = 24) was computed by the permutation test in imcRtools. Cyan bar represents P value < 0.05 whereas orange bars are not significant. I, Bar plot showing the change in spatial proximity of immune cell populations to NLRC5-high vs. NLRC5-low epithelial cells, averaged across all ROIs. J, Bar plot showing the change in spatial proximity of CD8 T cells to NLRC5-high vs. NLRC5-low epithelial populations, averaged across all ROIs. P values were determined by the sccomp differential composition test (C), the two-sided unpaired Wilcoxon test (E and H), and the two-sided paired Wilcoxon test (I and J). *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001.

Supplementary Figure S9. Identification of broad cell types in IMC data.

Supplementary Figure S5. Flow cytometry experiment details.

Bronchial premalignant lesions (PML), precursors of lung squamous cell carcinoma, have distinct molecular subtypes. The proliferative subtype, enriched with bronchial dysplasia, had decreased expression of an antigen-processing/presentation gene coexpression module in progressive/persistent versus regressive PMLs, suggesting a functional impact of these genes on immune evasion. In this study, we performed miRNA sequencing, miRNA in situ hybridization, and spatial proteomics of bronchial biopsies from patients at high risk for lung cancer. An miRNA-gene network analysis identified hsa-miR-149-5p as a potential regulator of the antigen presentation gene module. Staining on adjacent biopsy tissue showed that hsa-miR-149-5p was predominantly expressed in the epithelium and upregulated in progressive/persistent proliferative lesions. Targets of this miRNA, the transcriptional coactivator of MHC-I gene expression, NLRC5, and the genes it regulates, were downregulated in these lesions. Decreased NLRC5 expression reduced both IFNγ-induced MHC-I surface expression and CD8 T-cell cytotoxicity in lung squamous cancer cells. In PMLs, basal cells with high levels of NLRC5 were in close spatial proximity to CD8 T cells, suggesting that these cells exhibit increased functional MHC-I gene expression in vivo. These findings indicate a functional role for hsa-miR-149-5p in PML progression/persistence and suggest this axis as a potential therapeutic target for PML immunomodulation.

Supplementary Figure S1. Correlation between GSVA scores of miRNAs and gene modules from the miRNA-gene network.

Supplementary Figure S7. Performance of the hsa-miR-149-5p spot detection classifier.

Supplementary Figure S12. Representative IMC images.

Immune and stroma cell composition and epithelial NLRC5 expression is associated with PML progression. A and B, Volcano plots showing the differential composition of progressive/persistent (n = 7) compared with regressive (n = 7) PMLs for (A) cell populations identified in epithelium and (B) cell populations identified in stroma regions. C, Line plot showing the estimated marginal mean of NLRC5 expression with SD in epithelial cells of progressive/persistent (dark green) and regressive (light green) PMLs, respectively, across the three histology regions. Although histology groups were modeled as categorical, the line plot is used to facilitate visualization of trends across histology grades. D, Bubble plot showing the differential expression of NLRC5 between progressive/persistent and regressive PMLs in each epithelial cell population across the three histology contrasts. Color represents the estimate, and dot size represents the log P value computed by a linear mixed-effects model. P values were determined by the sccomp differential composition test (A and B) and linear mixed-effects models (C and D). *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001.

miRNA–gene network analysis identifies miRNAs associated with the antigen presentation gene module downregulated in progressive proliferative subtype PMLs. A, miRNA–gene module network plot shows 321 miRNAs connected to nine gene modules, in which gray circles are miRNAs and colored circles are the gene modules. Edges reflect a miRNA that may have strong suppressive regulatory role over the predicted target genes within the connected gene module based on the two statistical test filters. B, Pearson correlation densities between the module-associated miRNAs connected to each gene module in the miRNA–gene module network and the corresponding gene module GSVA score. C, miRNA–gene network plot of miRNAs connected to the antigen presentation gene module [module 9, purple circle in (A)]. Edges show the significant negative correlations between the miRNAs (gray circles) and the predicted target genes (purple circles). D, Boxplots of residual expression levels of the nine miRNAs specifically connected to module 9 in proliferative subtype biopsies stratified by their outcome status (progressive/persistent samples in dark green and regressive samples in light green). *, P < 0.05, determined by a linear mixed-effects model.

Supplementary Figure S4. Expression of hsa-miR-149-5p was up-regulated and expression of NLRC5 was downregulated in LUSC tumor compared to adjacent benign tissue.

Supplementary Figure S3. Expression level of hsa-miR-149-5p was significantly negatively correlated with that of NLRC5 in lung-related datasets.