Animal and housing care
Male C57Bl/6 mice,15 weeks in age were purchased from Charles River laboratories. Mice were housed separately in ventilated cages with normal chow and water ad libitum, in the Stanford Veterinary Service Center. The room was maintained at 23 °C, 60% relative humidity, and with a 12 h light/12 h dark cycle. Mice were continuously monitored for any possible health issues and body weights were measured thrice weekly. All procedures were performed in compliance with the ARRIVE guidelines 2.0 and was approved by the Stanford University Institutional Animal Care and Use Committee. All experimental procedures with mice were performed in accordance with the relevant guidelines and regulations as reviewed and approved by the Stanford University Institutional Animal Care and Use Committee (approval number: APLAC-31912). The study was designed and performed in compliance with the ARRIVE guidelines 2.0.
Experimental design
Mice were randomly allocated to control (n = 10), LD (n = 12) and HD CSE groups (n = 12) (Fig. 8A). All CS exposures were via whole body administration. LD mice were exposed to CS for 1 h/day. HD mice were exposed to 4 h of CS per day in 1-1.5 h increments, with room air breaks. Earlier studies in the airways and specifically in the larynx, have shown that 1-4 h/day of CSE is sufficient to induce histopathological changes 12,13,77,78. LD and HD CS exposures were performed for 5 days/week, M-F for a total of 4 weeks as per the subacute inhalation toxicity guidelines regulated by the Organization for Economic Co-operation and Development (OECD) 79. Before 4 weeks of CSE, the mice in both groups were acclimatized to ~20 min of CS per day, for 5 days (M-F). The exposure chamber concentration (610 µg/L) and delivered dose (mg/body weight) per cigarette (5.4882 mg/kg) were similar to our previous study 13. Control mice were left exposed to HEPA filtered room air conditions.
Cigarette smoke generation
Mainstream CS from Kentucky 3R4F reference cigarettes was generated using inExpose smoking equipment (SCIREQ, Montreal, QB, Canada), consistent with our previous studies 12–14. Exposure conditions were per the Federal Trade Commission (FTC)/International Standard Organization (ISO) standards (ISO 1991). Essential parameters were a puff duration of 2 s, a puff volume of 35 mL, a puff count of 9 per cigarette, and a puff frequency of 1 puff/min. Based on these puff characteristics, LD and HD group mice were exposed to 9 and 27 cigarettes/day respectively for 4 weeks. CS exposures were monitored in real-time using a probe, Microdust pro (Casella, UK) (Fig. 8B).
Animal euthanasia, tissue collection, and processing
Mice were euthanized by cervical dislocation under deep anesthesia by isoflurane 24 h after final CSE. Larynges were harvested and split between histological and RNA-related analyses. For histological analyses, larynges were fixed in 4% paraformaldehyde overnight at 4 °C and transferred to 70% ethanol (n = 5 control, n = 6 each in both CSE groups). HistoWiz, Inc. (http://www.histowiz.com; Brooklyn, NY, USA) processed tissues for paraffin embedding and performed coronal section microtomy per established methods 12–14,80,81. 5 µm coronal sections were used for immunohistochemical assessment of immune cell infiltrates. The remaining larynges across experimental groups (n = 5 control, n = 6 each in both CSE groups) were assigned for RNA extraction and stored in RNAlater solution (Catalog no. AM7020, Invitrogen, CA, USA) at -80 °C until RNA isolation.
Total RNA isolation and reverse transcription
Whole larynges were thawed and hemisected. Mucosal and submucosal tissues, starting from the glottis until the lower end of subglottic regions were harvested from the hemisected larynges. Total RNA was extracted according to the RNeasy Mini Kit (Qiagen, Germantown, MD, USA). RNA quantity and quality (A260/A280) were estimated using the Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Rockford, IL, USA). Samples from control (n = 1), low and high dose groups (n = 2 each) with poor RNA yield and A260/A280 values were eliminated from subsequent analysis. RNA integrity number (RIN) values were computed for remnant RNA elutes (n = 4 per group). Isolated RNA from all samples had an acceptable RIN value > 8 (Supplementary Fig. S4). Subsequently, reverse transcription was performed with 0.6 μg of RNA per reaction according to the manufacturer’s instructions in the iScript™ complementary DNA (cDNA) synthesis kit (Catalog no. 170-8890; Bio-Rad, Hercules, CA, USA).
Inflammatory gene expression assay and analysis with NanoString technology
Laryngeal mRNA was analyzed with a NanoString nCounter mouse inflammation v2 panel (Catalog no. XT-CSO-MIN2-12, NanoString Technologies, Seattle, WA, USA). Isolated RNA (~250 ng) from all biological samples across control, LD, and HD groups (n = 4 each) were specific to this panel. The number of occurrences of these unique bar-coded probes was counted by NanoString’s nCounter instrument for every sample. Raw data counts were obtained for further analysis.
NanoString’s data analysis application software, nSolver version 4.0 was utilized for the analysis of raw data counts. All experimental samples were annotated based on their exposure conditions. The background threshold was calculated as 12 using the recommended formula: Mean negative probe counts + 2 [standard deviation of negative probe counts]. Probes that had expression < 12 were eliminated from the analysis. The dataset was normalized using positive control and housekeeping probes. The resulting normalized and annotated grouped data were used for subsequent analysis. Global gene expression across different experimental conditions was assessed via complete linkage hierarchical clustering with Euclidean distance metrics using annotated grouped data in R. Differential gene expression was performed on the normalized dataset via advanced analysis module version 2.0 on nSolver. Benjamini-Hochberg false discovery rate (FDR) correction was employed to compute adjusted p-values. Comparisons between two independent groups a) HD vs control and b) LD vs control were conducted. Genes that passed a modest statistical cutoff 82 of p-value ≤ 0.05, FDR < 0.2, and fold change > +1.25 or fold change < -1.25 were identified to be significantly upregulated or downregulated.
PPI network analysis using StringDB and Cytoscape
The interaction of DEGs at the protein level was analyzed by using StringDB (string-db.org). Medium confidence PPI clustered networks with a minimum score of 0.4 and Markov Cluster Algorithm inflation parameter as 2 were built. KEGG pathway enrichment analysis was performed, with FDR ≤ 0.05 and minimum gene count as 3 considered to be statistically significant 83. Within the PPIs, the most interactive gene nodes, known as hub genes, were identified using the cytoHubba plugin 84 in the open-source software application, Cytoscape (http://www.cytoscape.org/; version 3.8.2). Out of the 12 different topological algorithms present in the cytoHubba plugin in the Cytoscape platform 84, MCC algorithm has a better execution in hub gene identification 85 and we used this topological algorithm to identify the top 10 hub genes.
qPCR validation of hub gene mRNA expression
qPCR was performed with CFX96 Real-Time PCR Machine (Bio-Rad) to validate NanoString gene expression levels of the key hub genes identified from PPI analysis. Mouse-specific primer sequences were designed based on NanoString target probe sequences using ApE, a plasmid editor (Table 1) 86. qPCR was carried out using iTaq™ Universal SYBR® Green Supermix (Catalog no.1725120; Bio-Rad, Hercules, CA, USA) in addition to template cDNA and specific primers. All qPCR experiments were conducted in duplicates. Relative quantitative analysis of gene expression was completed using the standard comparative Ct method (2-ΔΔCt) by using the recorded raw Ct values of the gene of interest and housekeeping gene 87. Using RefFinder, Gusb was selected as the most stable reference gene for qPCR normalization based on expression stability assessment among the other common candidate reference genes, Gapdh, Hprt, and Tubb5, present in the NanoString panel 88,89.
Upstream regulators and biofunctions identification via IPA
Upstream analysis and diseases and biofunctions modules were run on the DEGs via core analysis in IPA. All genes in the normalized NanoString dataset were used as the background reference set to compute enrichment. Upstream regulators and biological functions with Fisher’s exact p-valueof overlap < 0.05 and z-score >2 or z-score < -2 were identified to be significant.
Immunohistochemical assessment of immune cell infiltrates
Standard IHC procedures to stain F4/80+ macrophages, Myeloperoxidase/MPO+ neutrophils, and CD3+ T cells were performed by HistoWiz, Inc (histowiz.com). Stained slides were scanned by Histowiz using Aperio AT2 microscope (Leica Biosystems, Germany). These scanned virtual slides were used to perform an automated digital analysis of 3,3'-Diaminobenzidine (DAB) positive immune cell counts via QuPath version 0.3.0. For this analysis, mid membranous regions in the right and left vocal folds were first traced at 20x (Supplementary Fig. S5). Positive immune cell nuclei were counted in these regions by using the positive cell detection tool. In terms of the right and left subglottic regions, the entire subglottic regions (including epithelial and subglandular regions), starting from the lower end of VF until the first tracheal ring were traced at 10x (Supplementary Fig. S5). Positive immune cell nuclei were counted in these regions by using the positive cell detection analysis tool. Positive cells were detected using average brown DAB staining within the nucleus (Nucleus: DAB OD mean). The resulting data were expressed as % positive (number of positive detections / total number of cell detections in the traced area in µm2) averaged across the right and left sides of both VF and subglottis.
Statistical analysis
Statistical tests for body weight, IHC, and qPCR analysis were performed using GraphPad Prism version 9.0. Data comparisons were made between control, LD, and HD exposure groups by using a) two-way analysis of variance (ANOVA) for bodyweight pre/post CSE comparisons and b) unpaired two-tailed Student t-test for qPCR and IHC analyses. All bar graphs are represented as means ± standard error of mean (SEM). p-value ≤ 0.05 was considered statistically significant. All statistical analyses pertaining to NanoString global and differential gene expression, network analysis, and IPA are described above. Heatmap, volcano, dot and scatter plots were plotted using RStudio version 1.4.1717. Inter-platform correlation between NanoString and qPCR gene expression data of the selected key hub genes were performed via Pearson correlation in RStudio version 1.4.1717. The strength of correlation as indicated by Pearson correlation coefficient R was interpreted according to previously published guidelines 90.