Transcriptional profiling of lung macrophages following ozone exposure in mice identifies signaling pathways regulating immunometabolic activation

Macrophages play a key role in ozone-induced lung injury by regulating both the initiation and resolution of inflammation. These distinct activities are mediated by pro-inflammatory and anti-inflammatory/pro-resolution macrophages which sequentially accumulate in injured tissues. Macrophage activation is dependent, in part, on intracellular metabolism. Herein, we used RNA-sequencing (seq) to identify signaling pathways regulating macrophage immunometabolic activity following exposure of mice to ozone (0.8 ppm, 3 hr) or air control. Analysis of lung macrophages using an Agilent Seahorse showed that inhalation of ozone increased macrophage glycolytic activity and oxidative phosphorylation at 24 and 72 hr post exposure. An increase in the percentage of macrophages in the S phase of the cell cycle was observed 24 hr post ozone. RNA-seq revealed significant enrichment of pathways involved in innate immune signaling and cytokine production among differentially expressed genes at both 24 and 72 hr after ozone, while pathways involved in cell cycle regulation were upregulated at 24 hr and intracellular metabolism at 72 hr. An interaction network analysis identified tumor suppressor 53 (TP53), E2F family of transcription factors (E2Fs), Cyclin Dependent Kinase Inhibitor 1A (CDKN1a/p21), and Cyclin D1 (CCND1) as upstream regulators of cell cycle pathways at 24 hr and TP53, nuclear receptor subfamily 4 group a member 1 (NR4A1/Nur77), and estrogen receptor alpha (ESR1/ERα) as central upstream regulators of mitochondrial respiration pathways at 72 hr. These results highlight the complex interaction between cell cycle, intracellular metabolism, and macrophage activation which may be important in the initiation and resolution of inflammation following ozone exposure.

Methods
Total RNA was extracted as described above from 3 mice/treatment group. In a pilot study, we found that 3 mice were sufficient to identify a significant difference in Ptgs2 gene expression by qPCR at α = 0.05 and power = 80%. RNA integrity numbers (RINs) were confirmed to be ≥ 8.8 using a 2100 Bioanalyzer Instrument (Agilent, Santa Clara, CA). cDNA libraries were prepared using mouse TruSeq® Stranded Total RNA Library Prep kit (illumina, San Diego, CA) and quantified using a KAPA Library Quantification kit (Roche, Pleasanton, CA). cDNA libraries were sequenced (75 bp single-ended, ~35-44M reads per sample) on an Illumina NextSeq instrument. Raw reads in FastQ files were trimmed using Trimmomatic-0.39 (Bolger et al. 2014) and quality control of trimmed files performed using FastQC. Salmon was used to align reads in mapping-based mode with selective alignment against a decoy-aware transcriptome generated from mouse transcriptome GENCODE Release M23 (GRCm38.p6). Estimated counts per transcript were generated using the gcBias flag and normalized to transcript length to correct for potential changes in gene length across samples from differential isoform usage (Love et al. 2016; Patro et al. 2017). Transcript level quantitation data were aggregated to the gene-level using tximport (Soneson et al. 2015). Differential gene expression analysis was performed with air exposed mice as controls using DESeq2 with corrections for differences in library size (Love et al. 2014) in R version 4.0.3. Significantly enriched canonical pathways and upstream regulators were identified with Ingenuity IPA Version 65367011 (QIAGEN Inc, https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/) using a right-tailed Fisher’s Exact Test (Krämer et al. 2014). A less stringent criteria (fold change > 1.3 and experimental false discovery rate [padj] < 0.05) was used to augment the number of genes included in the pathway analysis (Bennett et al. 2024). Data were deposited NCBI’s Gene Expression Omnibus (Edgar et al. 2002) and are accessible through GEO Series accession number GSE237594 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE237594).

巨噬细胞（Macrophages）通过调控炎症的启动与消退，在臭氧诱导的肺损伤中发挥关键作用。这两类截然不同的功能由促炎型与抗炎/促消退型巨噬细胞介导，它们会依次在损伤组织中募集。巨噬细胞的活化在一定程度上依赖于细胞内代谢。本研究通过RNA测序（RNA-sequencing, RNA-seq），探究了小鼠暴露于0.8 ppm臭氧3小时或空气对照后，调控巨噬细胞免疫代谢活性的信号通路。

使用安捷伦Seahorse分析仪对肺巨噬细胞进行分析后发现，臭氧吸入可在暴露后24小时与72小时增强巨噬细胞的糖酵解活性与氧化磷酸化水平。臭氧暴露24小时后，可观察到处于细胞周期S期的巨噬细胞占比升高。RNA-seq结果显示，在臭氧暴露后24小时与72小时的差异表达基因中，先天免疫信号与细胞因子生成相关通路均显著富集；而细胞周期调控通路仅在24小时上调，细胞内代谢通路则在72小时上调。相互作用网络分析鉴定出，肿瘤蛋白53（tumor suppressor 53, TP53）、E2F转录因子家族（E2F family of transcription factors, E2Fs）、细胞周期蛋白依赖性激酶抑制剂1A（Cyclin Dependent Kinase Inhibitor 1A, CDKN1a/p21）与细胞周期蛋白D1（Cyclin D1, CCND1）是24小时细胞周期通路的上游调控因子；而TP53、核受体亚家族4A组1成员1（nuclear receptor subfamily 4 group a member 1, NR4A1/Nur77）与雌激素受体α（estrogen receptor alpha, ESR1/ERα）则是72小时线粒体呼吸通路的核心上游调控因子。本研究结果揭示了细胞周期、细胞内代谢与巨噬细胞活化之间的复杂相互作用，这或许在臭氧暴露后的炎症启动与消退过程中发挥重要作用。

## 实验方法
总RNA提取参照前述方法，每组取3只小鼠的样本。预实验结果显示，当显著性水平α=0.05、检验效能power=80%时，3只小鼠即可通过实时定量PCR（qPCR）检测到Ptgs2基因表达的显著差异。使用2100生物分析仪（Agilent, 圣克拉拉, 加利福尼亚州）确认RNA完整性指数（RINs）≥8.8。使用小鼠TruSeq® 链特异性总RNA文库制备试剂盒（Illumina, 圣迭戈, 加利福尼亚州）构建cDNA文库，并采用KAPA文库定量试剂盒（Roche, 普莱森顿, 加利福尼亚州）进行文库定量。随后在Illumina NextSeq测序仪上进行测序（75 bp单端测序，每个样本约3500万至4400万条reads）。

使用Trimmomatic-0.39对FastQ格式原始reads进行质控修剪（Bolger等, 2014），并通过FastQC对修剪后的reads进行质量评估。采用Salmon软件以映射模式结合选择性比对策略，将reads比对到由小鼠转录组GENCODE Release M23 (GRCm38.p6)构建的含诱饵序列的转录组中。使用gcBias参数生成转录本的估计计数，并针对样本间因差异异构体使用导致的基因长度变化进行转录本长度归一化（Love等, 2016; Patro等, 2017）。使用tximport工具将转录本水平的定量数据汇总为基因水平数据（Soneson等, 2015）。以空气暴露小鼠为对照，使用R版本4.0.3中的DESeq2软件进行差异基因表达分析，并对文库大小差异进行校正（Love等, 2014）。使用Ingenuity IPA Version 65367011（QIAGEN公司, https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/）进行右侧单侧Fisher精确检验，以鉴定显著富集的经典通路与上游调控因子（Krämer等, 2014）。为扩充通路分析纳入的基因数量，采用较为宽松的筛选标准（折叠变化>1.3，实验校正后假阳性率[padj]<0.05）（Bennett等, 2024）。

本研究数据已提交至NCBI基因表达综合数据库（Gene Expression Omnibus, GEO）（Edgar等, 2002），可通过GEO系列登录号GSE237594（https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE237594）获取。