依赖于NRF2 的lncRNA-mRNA 共表达网络的全基因组识别及其在非小细胞肺癌进展中的作用

唐建军，韩春宾，纪玉龙，张剑锋，姜斯聪

（1.南昌大学第一附属医院呼吸与危重症医学科，江西 南昌 330052；2.江西省肿瘤医院胸外科，江西 南昌 330029；3.江西省肿瘤医院转化肿瘤研究重点实验室，江西 南昌 330029；4.江西省肿瘤医院病理科，江西 南昌 330029）

在全球所有癌症类型中，肺癌（lung cancer）的死亡率最高[1]。因此，了解关键致癌驱动因素对疾病生物学和临床结果的影响对于设计新疗法至关重要。非小细胞肺癌（non-small cell lung cancer，NSCLC）和小细胞肺癌是肺癌的亚型[2]。NSCLC 被细分为离散的组织学和基因组亚型，每个亚型都表现出独特的生物学和临床特征[1]。Kelch 样环氧氯丙烷相关蛋白-1（KEAP1）和核因子E2 相关因子2（NRF2）基因位点的突变是NSCLC 中的常见事件，并导致转录因子NRF2（抗氧反应的关键调控因子）的构成稳定[3，4]。NRF2 是碱性亮氨酸拉链结构域蛋白家族的成员，可调节抗氧化酶和含有抗氧化反应元件的细胞保护蛋白基因的表达[5，6]。活性氧（ROS）是氧代谢的部分还原氧副产物[7，8]。ROS 的积累可导致DNA 氧化损伤并促进表观遗传改变，从而导致驱动基因失调和癌症发病机制[9，10]。NRF2 在对损伤和炎症的抗氧化反应中减轻氧化应激起关键作用[10]。它受负调节因子KEAP1 蛋白的严格调控，并可与NRF2 蛋白结合，诱导其泛素化和降解[11]。然而，关于NRF2 和KEAP1 在癌症发病机制中的作用还缺乏共识。NRF2 被认为是肺癌的致癌基因且与NSCLC的发病机制和进展有关。并且，肺癌组织中的NRF2上调也与不良治疗预后相关[11-14]。此外，肺癌的诱发机制似乎不同于在其他器官系统和肺部疾病中观察到的机制，如肺气肿[15]、高氧[16]和呼吸道合胞体病毒[17]，这些都是KEAP1 依赖性的。长链非编码RNA（lncRNA）长度＜200 个核苷酸，调控基因表达和关键的生物过程[18]，包括细胞周期、细胞凋亡[19-21]以及染色质修饰和重塑。鉴于它们影响着重要的生物过程和许多疾病，增加对其遗传学和生物学的了解已成为近期研究的新热点。研究表明[23，24]，lncRNAs 在肺癌发病机制中调节NRF2 表达中的作用。为此，本研究主要探究NRF2 调控的lncRNA 在功能上与肺癌有关信号通路的关系，揭示在肺癌中NRF2 对lncRNA 调控的机制。

1 资料与方法

1.1 资料来源 在TCGA数据库（https://can cergenome.nih.gov/）获取486 个NSCLC 样本和50个正常样本的RNA-seq 数据。使用FastQC（http://www.bioinformatics.babraham.ac.uk/projects/fastqc）检查原始数据质量，通过FASTX-Toolkit v.0.0.13（http://hannonlab.cshl.edu/fastx_toolkit/）去除低质量碱基。从GEO 数据库（GSE118844）中提取公共NRF2 依赖性转录组和H3K27ac 染色质免疫沉淀测序数据，H3K27ac ChIP-seq 数据包括3 组NRF2-KD 肺癌A549 细胞系和3 组用于免疫沉淀的对照细胞，以及两组DNA 输入（基因组DNA 作为阴性对照）。RNA-seq 数据来自KEAP1 突变条件（NRF2 激活）下的NRF2-KD 和A549 细胞。

1.2 染色质免疫沉淀测序（ChIP-seq）使用Bowtie 2软件分析公开可用的ChIP-seq 数据，将全基因组ChIP-seq 读数与人类GRCh38 参考基因组对齐。错误发现率设置为q值＜0.05（Benjamini-Hochberg 调整的P值）。ChIP-seq 算法（MACS v.14 软件）的模型分析被用于识别NRF2 依赖的富含H3K27ac 的ChIP-seq 峰；控制样本被设置为输入。接下来，使用“bedtools cluster”对每个样本中的ChIP-seq 峰进行聚类。在TSS 上游10 kb 和下游3 kb 的区域中发现了ChIP-seq 启动子相关的峰相关基因。

1.3 转录组测序 使用TopHat2[25]将高质量数据与GRCH38 参考基因组进行对比，并最多允许4 个错配（Ensembl 释放61 预测基因模型）。使用外显子级估量来改进差异基因表达，通过计算映射到基因外显子的读数量来获得每个基因的计数。通过edgeR 软件[26]使用以下标准来识别DE lncRNA：倍数变化（FC）≥1.5 或FC≤0.67；和错误发现率（FDR，P值≤0.05）。

1.4 LncRNA 预测和方向识别 为了系统地分析lncRNA 表达模式，采用管道法鉴定lncRNA[27]。该程序是根据Cufflinks 软件设计的[28]。

1.5 功能富集分析 使用KOBAS2.0 服务器识别GO项、KEGG 和Reactome 路径[29]。采用超几何检验和Benjamini-Hochberg FDR 调整的P值来定义各术语的富集度。

1.6 LncRNA-mRNA 共表达网络分析 使用CNC 网络来计算TCGA NSCLC 样本（n=550）中共有lncRNAmRNA 对的PCC（临界值≥0.4）。使用Cytoscape3.7.2版软件（The CytoscapeConsortium，SanDiego，CA，USA）对lncRNA-mRNA 网络进行可视化。

1.7 Kaplan-Meier 生存分析 对不同的上调和下调组进行Kaplan-Meier 生存分析并通过R 包（survival and survminer）（https://cran.r-project.org/web/packages/survival/index.html）生存分析代码进行可视化，组间比较采用数秩检验（Log-rank），TCGA NSCLC 数据集中的中位基因表达设置为上调和下调组的截止值。

1.8 统计学方法 在使用每百万份转录本标准化方法后，将样本中每个基因的读取（内部脚本）用于下一代测序数据的可视化和基因组注释。使用R 3.6.1（https://www.r-project.org/）进行统计学分析，两组间比较采用独立样本t检验。

2 结果

2.1 肺癌中NRF2 的过度表达和 KEAP1 突变体A549 细胞中NRF2 结合靶点的鉴定 从486 例癌症和50 例正常样本的TCGA 数据集中发现NRF2 在肺癌患者中高表达（图1A）。NCBI/GEO/H3K27ac/ChIP-seq（GSE118840）数据集包含3 个含有NRF2 siRNA（NRF2-KD）的NRF2 敲低A549 样本，3 个用于免疫沉淀的含有对照siRNA 的对照样本和两个以对照siRNA 作为输入的样本（基因组DNA）。与NRF2-KD 样品相比，从3 个对照样品中特异性鉴定出8277、3808 和6858 个富集峰；而与对照样品相比，在3 个NRF2-KD 样品中共发现2645、5396和4020 个峰（图1B）。H3K27ac 中超过80%的结合峰的宽度小于2151 bp（图1C）。对NRF2NRF2 依赖的H3K27ac 峰进行文氏图分析，并从3 个ChIP-seq重复中至少2 个中确定了3006 个候选NRF2 依赖的H3K27ac 峰（图1D）。KEGG 富集显示，NRF2 具有发挥多方面的潜在作用（图1E）。依赖于NRF2 的H3K27ac 沉积峰在转录起始位点（TSS）附近的10 kb 区域显示出较高的比例（15%）（图1F）。峰的全基因组分布显示114 个lncRNA、10 个sense_intronic、124 个假基因、3 个sense_overlapping、102 个反义、34 个microRNA（miRNA）、28 个misc_RNA、452 个编码蛋白，以及71 种其他基因类型（图1G）。对NRF2 依赖的H3K27ac 峰的基因本体论（GO）分析表明，该峰富集了各种代谢过程（图1H）。对照样本和NRF2-KD 样本中有代表性的lncRNAs 基因组快照（RP11-445P17.8）显示，与NRF2-KD 样本相比，对照样本具有强的H3K27ac 结合峰（图1I），表明肺癌中上调的NRF2 可以在lncRNA 位点形成H3K27ac 标记的增强子。

图1 KEAP1 敲低A549 细胞中NRF2 结合靶点的识别

2.2 NRF2-KD 与具有KEAP1 突变的对照A549 细胞之间的DE lncRNA 分析 结果表明，对照组在TSS 附近有很强的结合模式，而NRF2-KD 组在TSS上的结合微乎其微（图2A、2B，图3A、3B）。火山图和热图显示，NRF2-KD 与对照样品中有368 个lncRNA 基因表达下调，404 个表达上调（图3C、3D）。使用550 个LUSC 样本的TCGA 数据集通过共表达研究lncRNA 和mRNA 表达之间的关系，DE lncRNA 和共表达mRNA 的GO 和KEGG 功能富集分析见图3E、3F，表明DE lncRNA 和mRNA 在肺癌进展中发挥作用。

图2 NRF2-KD 与具有KEAP1 突变的对照A549 细胞之间的DE lncRNA 分析

图3 结合DE mRNA 功能途径的分析

2.3 NRF2 介导的lncRNA 在肺癌中的调控作用 为了建立NRF2 调控的lncRNA 与肺癌之间的临床关系，从数据库中筛选已知的lncRNA，并预测罕见的新lncRNA。GO 和KEGG 分析确定依赖于NRF2 的DE mRNA 和H3K27ac 结合的mRNA 的重叠基因富集途径，共发现一组由NRF2 依赖性增强子调节的18 个lncRNA（图4A、4B）。这18 个lncRNA 表达谱的热图显示了两个不同的簇，一个对应于对照组，另一个对应于NRF2-KD 组。然后，对TCGA LUS 数据集中的18 个lncRNA 进行共表达分析，识别它们相关的mRNA。共表达mRNA 分析见图4C、4D，DE lncRNA-mRNA 通路图见图4E。

图4 NRF2 介导的LncRNA 在肺癌中的调控作用

2.4 lncRNA LINC00488 对肺癌患者预后的影响 对18 个NRF2 依赖性H3K27ac 结合的lncRNA 进行Kaplan-Meier 生存分析，结果显示18 个lncRNA 中LINC00488 上调与预后不良相关。ChIP-seq 结果和RNA-seq 结果中LINC00488 的读取分布见图5A、5B、5C；ChIP -seq 分析识别出lncRNA RP11.227H15.4（Chr10:69231262 -69233332）的NRF2 依赖性增强子区域（图6A）。此外，Kaplan-Meier 生存曲线显示，RP11.227H15.4 基因表达上调的LUSC 患者预后不良（图6B）。

图5 lncRNA LINC00488 对肺癌患者预后的影响

图6 调节lncRNA（RP11.227H15.4）对肺癌患者预后的影响

3 讨论

肺癌是一种异质性疾病，了解关键致癌驱动因素对其生物学和患者临床结局的影响对于设计新疗法至关重要。NSCLC 是最普遍的肺癌类型[30]，lncRNA 参与肿瘤发生是一个新兴的主题[31，32]。LncRNA 可能通过改变基因表达发挥与经典癌基因或肿瘤抑制基因相似的作用[33]。NRF2 在调节细胞内氧化微环境、炎症反应和细胞稳态中起关键作用[34]。此外，它还可以在癌症中发挥致癌或保护作用[35-37]。

小的非编码RNA 在增殖、干细胞自我更新、细胞凋亡和化学/放射抗性等生理过程中发挥重要作用[38，39]。这些非编码RNA 也被证明可以在mRNA 和蛋白质水平上调节癌基因和肿瘤抑制基因。类似地，NRF2 可以通过miRNA 和lncRNA 在遗传和表观遗传水平上对癌基因和肿瘤抑制基因进行调节[40，41]。许多研究报道了NRF2 和lncRNAs 在肿瘤进展中的关系。如前列腺癌进展中的NRF2 表达受lncRNA TUG1 调节[42]。此外，吸烟和癌症相关的lncRNA-1/肺癌相关转录物1（SCAL1/LUCAT1）的表达受NRF2 在烟草烟雾诱导的肺癌中的调节[43]。此外，敲低NRF2 可以降低lncRNA 核缺氧调节的NLUCAT1 的表达；而NRF2 的过表达可以提高NLUCAT1 的表达并促进肺腺癌中的顺铂耐药[44]。因此，探究NRF2-lncRNA 的相互作用至关重要。尽管lncRNA 与调节NRF2 基因表达有关，但目前仍缺乏关于其在肺癌中作用的详细全基因组视图[45]。

本研究确定了一个依赖DE NRF2 的lncRNA网络，并通过全基因组RNA-seq 分析以及GO 和KEGG 通路富集分析探索它们在NSCLC 中可能的临床意义，结果显示其主要涉及的通路有：细胞色素P450 介导的异生物代谢，硫酸角质素糖胺聚糖生物合成和细胞色素P450 介导的药物代谢、代谢途径，丙氨酸、天冬氨酸和谷氨酸代谢，淀粉和蔗糖代谢，类固醇激素生物合成，化学致癌，卟啉和叶绿素代谢。值得注意的是，对A549 细胞（GSE118844）中NRF2 依赖性H3K27ac ChIP-seq 数据和NRF2-KD RNA-seq 数据的系统分析发现了368 个lncRNA 基因表达下调，404 个表达上调，并且共表达的mRNA具有自身免疫性疾病介导的细胞毒性、自噬调节、IgA 产生的肠道免疫网络、RIG-Ⅰ样受体信号通路和细胞粘附分子的富集。至少2 个DNA 样本和3006 个候选NRF2 依赖的H3K27ac 峰。该峰富集了各种代谢过程，如小分子代谢、受体介导的内吞作用、硫酸角质素生物合成、有丝分裂的G1/S 转换、水溶性维生素代谢和小GTPase 介导的信号转导的调节。本研究表明，NRF2 上调可以诱导细胞周期从G0/G1期向S 期转变，促进人支气管上皮细胞的肿瘤发生[46]，以及NRF2 在抗氧化防御、NADPH 的再生和其他代谢途径中的作用[13，23]。本研究发现了18 个lncRNA 共表达mRNA，表明这些基因可能与RNA剪接、DNA 依赖性DNA 修复、mRNA 加工、转录调控、有丝分裂、DNA 复制和转录激活可能关联。对共表达mRNAs 的KEGG 分析预测了与剪接体、谷胱甘肽代谢、赖氨酸降解、细胞色素P450 介导的异生物质代谢、同源重组、细胞周期、RNA 转运、泛素介导的蛋白水解和胞质铁硫簇组装途径的关联。最后发现lncRNA LINC00488 的高表达与肺癌患者预后不良有关。而lncRNA LINC00488 与许多癌症有关，如通过miR-376a-3p/PON2 的甲状腺癌[47]和通过miRNA-485-5p 的食道癌[48]。未来的研究将探索利用功能性敲入和敲除实验阐明NRF2-lncRNA LINC00488 在肺癌中的调控机制。

综上所述，本研究发现lncRNA LINC00488 可能在肺癌中受到NRF2 的潜在调节，其可能作为新的预后生物标志物和肺癌治疗靶点的潜在作用。但需要进一步的实验来阐明lncRNA 对NRF2 表达及其下游靶点在癌症进展中的互补作用。