肝豆状核变性(hepatolenticular degeneration)，又名威尔森疾病(Wilson disease, WD)，是一种常染色体隐性遗传性铜代谢障碍疾病。铜离子转运ATP酶β肽(ATPase Cu²⁺ transporting beta polypeptide，ATP7B) 基因突变是肝豆状核变性的主要致病原因^[1-2]。典型的肝豆状核变性通常伴随肝脏受累的表现(包括转氨酶升高、急/慢性肝炎、肝硬化和肝功能衰竭) 或神经系统症状(如运动障碍、肌张力障碍和精神状态等)^[3]。然而，部分患者在疾病发生的早期并不具备典型的肝豆状核变性症状，甚至约20%的患者检测不到ATP7B基因的纯合或复合杂合突变^[4-5]。这为疾病的早期诊断带来困难。

目前，国内外关于肝豆状核变性致病机制的基础研究多局限于ATP7B基因突变谱和突变热点的发现，较少考虑转录调控对疾病发生发展的影响。然而，越来越多的证据表明，转录调控和表观遗传机制在研究肝豆状核变性表型发生和发展中具有重要意义^[6]。到目前为止，仅有极少数的研究报道使用小鼠模型对肝豆状核变性的调控模式进行了初步的探索，研究人员利用小鼠肝脏组织分析了差异lncRNA和mRNA的表达，并讨论了转录调控和信号通路的变化^[7-8]。此外，有研究发现，ATP7B基因缺失的肝细胞通过激活细胞自噬降低细胞凋亡水平，能够缓解过量铜诱导的肝细胞损伤^[9]。然而，在疾病发生的过程中，与肝豆状核变性密切相关的转录调控因子和疾病本身涉及的转录调控机制尚不清楚。本研究基于转录组测序数据和生物信息学分析方法探讨了铜诱导的ATP7B基因敲除HepG2细胞系转录调控的分子机制，筛选出肝豆状核变性转录调控网络中的关键转录因子，为揭示疾病潜在的调控模式和促进疾病的分子功能机制研究提供了参考依据。

以“Wilson disease”为关键词检索GEO数据库(https://www.ncbi.nlm.nih.gov/geo/)，根据以下几条具体原则进行数据筛选工作：(1) 数据来源为人源；(2) 基于高通量转录组测序；(3) 至少3个生物学重复。最终筛选到GSE107323数据集进行后续分析^[9]。该数据集共包括12例样本，包括野生型HepG2细胞系和经ATP7B基因敲除的HepG2细胞系，以及铜诱导的野生型HepG2细胞系和铜诱导的ATP7B基因敲除的HepG2细胞系各3例，其测序平台为GPL18573。

分别利用R语言(version 4.0.3) 的limma^[10]、edgeR^[11]和DESeq2^[12] 3个软件包对上述野生型和经ATP7B基因敲除的HepG2细胞系进行差异表达基因(differentially expressed genes, DEGs) 分析，设置同时满足P < 0.05、|log₂FC| > 1且FDR < 0.05为DEGs的筛选阈值。对3个软件包的分析结果取交集，获得高可信度的、经基因敲除产生的DEGs。再采用同样的计算方法和阈值设定，对铜诱导的野生型Hep2G细胞系和铜诱导的ATP7B基因敲除的HepG2细胞系进行DEGs分析，获得可信度高的且铜诱导的DEGs。对两类差异基因取补集，扣除ATP7B基因敲除导致的差异基因，最终获得因铜诱导产生的差异表达基因。并使用R语言pheatmap软件包(version 1.0.12. https://CRAN.R-project.org/package=pheatmap) 绘制差异基因热图。

利用基因本体论(gene ontology, GO) 数据库对上述DEGs进行GO功能富集分析，包括生物学过程(biological process, BP)、细胞组分(cellular component, CC) 和分子功能(molecular function, MF) 3个方面。京都基因与基因组百科全书(Kyoto encyclopedia of genes and genomes, KEGG) 通路富集分析用于识别DEGs参与的主要生物学通路。本研究使用R语言clusterProfiler^[13]软件包对上述筛选出的DEGs进行GO和KEGG富集分析，设定P < 0.05并利用Benjamini-Hochberg法进行P值校正，识别显著富集的GO term和KEGG通路。

本研究利用Barabási等^[14]构建的包含了15 970个蛋白共217 160种相互作用的蛋白-蛋白相互作用(protein-protein interaction, PPI) 网络(命名为Barabási-PPI)，分别提取上、下调DEGs和Barabási-PPI网络中度分布排名在前5%节点的交集，定义为上、下调关键基因。然后，提取上、下调关键基因以及它们在Barabási-PPI网络中的邻居节点，作为疾病的关键功能模块。

进一步对上述筛选出的关键功能模块中的基因进行GO和KEGG富集分析，研究模块中基因的功能特性，具体方法同1.3节中的研究方法，设置P < 0.05为阈值并采用Benjamini-Hochberg法进行P值校正，对富集结果进行筛选。

Cistrome DB数据库是目前最全面的研究ChIP-seq和DNase-seq的数据库，主要收录了人和小鼠的转录和表观遗传基因调控数据^[15]。从该数据库中下载了与HepG2细胞系相关的转录因子列表，将转录因子和DEGs取交集筛选出差异表达转录因子。并利用R语言Hmisc软件包(version 4.5-0. http://CRAN.R-project.org/package=Hmisc) 中的rcorr函数计算差异表达转录因子和差异表达基因之间的相关性，定义cor > 0.5且P < 0.001为正相关，cor < –0.5且P < 0.001为负相关。最后，将节点度分布大于5的差异表达转录因子和差异表达基因导入Cytoscape^[16] (version 3.6.1) 中绘制出肝豆状核变性转录调控网络图。

利用R语言limma、edgeR和DESeq2软件包，对野生型和经ATP7B基因敲除的HepG2细胞系，以及铜诱导的野生型和铜诱导的ATP7B基因敲除的HepG2细胞系分别进行差异表达基因分析。对3个软件包的结果取交集后，再对这两类DEGs取补集，共筛选出1 034个由铜诱导产生的DEGs (图 1A)，其中509个基因表达下调，525个基因表达上调，并进一步绘制出差异表达基因热图(图 1B)，可见铜诱导的野生型和铜诱导的ATP7B基因敲除型细胞系的mRNA表达出现显著差异。

图 1 DEGs筛选 Fig. 1 Screening of DEGs. (A) Venn diagram showing overlap of DEGs by limma, edgeR and DESeq2 packages. (B) Heatmap of DEGs between wild-type and ATP7B-knockout HepG2 cells exposed to copper. WT1+Cu~WT3+Cu: wild-type HepG2 cells with Cu exposure; ATP7B-ko1+Cu~ATP7B-ko3+Cu: ATP7B-knockout HepG2 cells with Cu exposure; blue scale: down-regulated genes; red scale: up-regulated genes.

图选项

为了对DEGs影响的生物学功能和通路进行描述，分别对上、下调DEGs进行GO和KEGG富集分析，结果如图 2A–D所示。下调DEGs主要富集在细胞周期、细胞自噬、DNA修复和线粒体凋亡等GO term上(图 2A–B)。而上调DEGs则显著富集到糖代谢、小分子代谢、脂质代谢、醇代谢等GO term上，以及碳代谢、糖代谢等、丙酸盐代谢信号通路上(图 2C–D)。值得注意的是，从总体水平上看，上、下调DEGs显著富集的GO-CC和KEGG信号通路较少，并且在GO分子功能方面没有显著富集到任何术语。因此，需要进一步对这些DEGs进行研究。

图 2 差异表达基因的GO和KEGG富集分析 Fig. 2 GO and KEGG enrichment analysis of DEGs. GO enrichment results of down-regulated genes (A) and up-regulated genes (C); KEGG enrichment results of down-regulated genes (B) and up-regulated genes (D).

图选项

根据方法1.4中的描述筛选出上、下调关键基因和关键功能模块，分别得到关键下调基因37个(表 1) 和关键上调基因17个(表 2)。上、下调关键功能模块分别包含了3 785个和3 931个基因。为了进一步研究关键功能模块中基因的生物学功能和通路富集情况，分别对这2个关键功能模块中的基因进行GO和KEGG富集分析，富集结果如图 3A–D所示。值得注意的是，在DEGs的富集结果中，仅仅是差异表达基因无法显著富集到更多相关的生物学功能和信号通路，但通过PPI网络分析后，关键功能模块中的基因均能够显著地富集到某些GO term和KEGG通路上。其中，GO富集分析结果显示，上、下调关键功能模块中的基因主要定位于细胞-基质连接、染色体区域、剪接复合体、核糖体亚基、胞质核糖体等区域(图 3A、3C)；在生物学过程方面，DEGs主要参与mRNA加工、组蛋白修饰、RNA剪切、DNA代谢的调节、胞质翻译、蛋白磷酸化、蛋白定位、核糖核蛋白复合物发生(图 3A、3C)；在分子功能方面，DEGs主要参与了蛋白激酶活性、DNA转录因子结合、转录共调控活性、泛素样蛋白连接酶结合、钙粘蛋白结合等GO term (图 3A、3C)。KEGG分析结果显示，DEGs主要富集在乙型肝炎、MAPK信号通路、病毒性致癌、细胞衰老和凋亡、神经营养信号通路、神经变性途径、FoxO信号通路、T细胞受体信号通路、线粒体自噬等相关信号通路(图 3B、3D)。由于上、下调关键功能模块中存在部分重叠的基因，因此关键功能模块富集分析得到的GO term和KEGG通路存在部分重叠。但无论从上调功能模块的角度还是下调功能模块的角度分析，这部分GO term和通路均受到了影响，表明这些GO术语和KEGG通路对于肝豆状核变性具有重要的生物学意义。

图 3 关键功能模块中基因的GO和KEGG富集分析 Fig. 3 GO and KEGG enrichment analysis of genes in key functional modules. GO enrichment results of genes in down-regulated key functional modules (A) and in up-regulated key functional modules (C); KEGG enrichment results of genes in down-regulated key functional modules (B) and in up-regulated key functional modules (D).

图选项

为了减少对单个基因过度关注的局限性，本研究从整体角度系统分析了肝豆状核变性差异表达基因的转录调控方式，详细讨论了差异表达转录因子和差异表达基因之间的调控关系，并绘制了转录调控网络图(图 4)。最终得到的网络共包括11个差异表达转录因子和96个差异表达基因之间的547个调控关系。其中，U2AF1、NFRKB、FUS、MAX、SRSF1、CEBPA和RXRA这7个转录因子位于调控网络的中心，这些核心差异表达转录因子共同调控了大多数DEGs的表达(表 3)。此外，图 4显示，肝豆状核变性关键基因ATP7B还同时受到了差异表达转录因子HMG20A、NFATC3、ZNF792和HNF4G的调节。

图 4 差异表达转录因子和差异表达基因调控网络 Fig. 4 Differentially expressed transcription factors and differentially expressed genes regulatory network. Green square represents differentially expressed transcription factors, and yellow circle represents differentially expressed genes. The size of a node represents its degree. Red lines represent up-regulated, and blue lines represent down-regulated.

图选项

随着系统生物学方法的发展以及相关实验数据的积累，基于转录调控网络识别复杂疾病的致病基因和基因与靶点之间的调控关系，可以为基础研究提供潜在的候选靶基因。肝豆状核变性的总体发病率约为1/30 000，在亚洲人群中更为高发^[17]。该病虽然不是常见疾病，但在亚裔人群中并不罕见。而且肝豆状核变性是为数不多的几种可以治疗的遗传病之一，因此早期诊断和治疗对于患者意义重大。本研究利用转录组测序技术和生物信息学分析方法，分析了野生型Hep2G细胞系和经ATP7B基因敲除的HepG2细胞系之间的差异表达基因，以及经铜诱导的野生型Hep2G细胞系和经铜诱导的ATP7B基因敲除的HepG2细胞系之间的差异表达基因。保留铜诱导产生的差异表达基因，并对这部分DEGs进行GO和KEGG富集分析。其次，以Barabási-PPI网络为背景，进一步分析了肝豆状核变性关键基因和关键功能模块。最后，基于差异的表达基因和转录因子构建了转录调控网络，探讨了核心差异表达转录因子对肝豆状核变性的影响以及相关表型产生的可能原因，具有一定的参考价值。

本研究对关键功能模块中的基因进行了GO和KEGG富集分析，结果显示上、下调关键功能模块共同富集在某些GO term和KEGG通路上，如乙型肝炎、MAPK信号通路、FoxO信号通路、T细胞受体信号通路和神经营养信号通路上，以及mRNA加工、剪接、组蛋白修饰、蛋白磷酸化和细胞衰老凋亡等与细胞活动和功能相关的GO term上。从富集的KEGG通路结果来看，WD和乙型肝炎的症状之间很可能存在某些联系。并且临床证据也表明，WD的肝脏症状主要包括慢性隐性肝炎和急性暴发性肝炎^[2]。此外，Brady等发现，用于治疗WD的铜螯合剂也可用于治疗含有BRAF突变的癌症，而BRAF激酶磷酸化会刺激MAPK信号通路^[18]。MAPK信号通路在肝脏炎症的凋亡中具有重要意义，并且该信号通路通过调节肝星状细胞的活化、增殖和凋亡影响肝脏纤维化的形成^[19]。多项临床研究发现，WD患者中出现大量以中枢神经系统症状为首发症状的患者^{[2-3, 20]}，而本研究中关键功能模块内的基因显著富集到了神经营养信号通路和神经变性途径等信号通路上，提示差异基因及其关联的互作用基因的整体变化或许影响了WD神经表型的产生。Polishchuk等^[9]探讨铜诱导对野生型和敲除型HepG2细胞系的影响，经DAVID在线工具富集分析发现与野生型细胞相比，铜诱导对ATP7B缺陷的细胞系有更强的自噬激活作用。本研究侧重于比较铜诱导后的野生型和ATP7B缺陷型细胞系之间的转录调控差异，而clusterProfiler包的富集结果多集中在mRNA加工、剪接和翻译等与转录调控相关的GO term上。

FUS是一个DNA/RNA结合蛋白，FUS基因突变会造成神经元细胞死亡进而引发神经退行性疾病的发生^[21]。多项研究表明FUS基因突变与肌萎缩侧索硬化症相关，该病的主要症状表现为肌无力、肌萎缩和震颤痉挛^[21-22]。临床中，存在很大一部分WD患者，以震颤和肌张力障碍等神经症状为最常见的神经学特征。据文献报道有多达55%的肝豆状核变性神经系统表型患者在确诊时有震颤症状^[23]。在本研究结果中，FUS基因作为核心差异表达转录因子之一，提示FUS基因和它调控的基因的异常表达可能与肝豆状核变性神经表型的产生存在某种关联。

U2AF1和SRSF1同属于U2核糖核蛋白体依赖型剪接体相关基因，主要参与前体mRNA的剪接、加工以及调控细胞周期和凋亡。其中，SRSF1基因被认为是一个潜在的原癌基因，通过选择性剪接能够促进肝癌的发展^[24]。此外，还有研究表明SRSF1不仅可以调控iNKT细胞的发育分化，还能够增强对急性肝损伤的耐受力，减弱肝损伤程度，在其功能方面具有至关重要的作用^[25]。另外，Sanders等的研究结果表明，转录因子MAX在肝脏的发育、增殖和再生过程中都具有重要的调控作用^[26]。与kappaB结合蛋白(nuclear factor kappa-B, NFRKB) 相关的核因子在正常肝细胞中低/中表达，而在肝癌细胞的细胞核中大量积累，NFRKB的高表达提示不良的预后和较低的生存率，该基因的功能障碍与肝癌发展密切相关^[27]。

正常肝脏受到损伤时能够完全再生，而衰老的肝脏其增殖反应则明显减弱。CEBPA作为一种肝特异性蛋白，研究表明它能够通过抑制细胞周期依赖性激酶从而抑制肝细胞的增殖。Iakova等^[28]通过实验证实，衰老能够关闭CEBPA抑制E2F转录的肝脏特异性途径，从而导致肝脏增殖反应的丧失。本研究中，转录因子CEBPA的异常很可能影响了WD患者肝脏的再生增殖能力。另一方面，肝星状细胞在肝纤维化中起着重要作用，CCAAT增强结合蛋白(CCAAT enhancer binding proteins, CEBPs)，特别是CEBPA控制了脂肪细胞的分化。研究表明，通过CEBPA转染激活的肝星状细胞后，能够同时上调RXRA的表达，表明CEBPA和RXRA与肝脏脂肪细胞的分化和肝纤维化存在密切关联^[29]。此外，RXRA能够与受体基因形成异二聚体，在配体的作用下，对靶基因的表达进行调控^[30]。WD患者的肝细胞会出现不同程度的脂肪变性^[31]，而Kim等的研究结果表明，RXRA等相关基因表达的改变影响了肝脂肪代谢的改变^[32]。并且，Meacham等以肝豆状核变性小鼠为研究对象，探讨了铜-锌平衡对小鼠肝脏的影响，其研究结果表明，以RXRA为代表的一系列核受体基因与WD的转录调控密切相关，这些基因共同参与了细胞的分化、发育，以及脂质、蛋白质代谢等细胞活动^[8]。Meacham等的研究结果与本研究关于RXRA作为核心转录因子调控疾病的研究结论基本一致。

ATP7B基因作为肝豆状核变性的关键基因，除了受到7个核心差异表达转录因子的调控外，HMG20A、NFATC3、ZNF792和HNF4G也共同调节了该基因的表达变化。其中，Wooton-Kee等^[33]报道了肝豆状核变性患者由于受到过量铜的特异性和选择性影响，导致以RXR、HNF4等为例的一系列核受体基因与启动子元件的结合减少，使核受体靶基因mRNA表达减少的现象，这与本研究的相关结果一致。

肝豆状核变性作为少数可以有效控制的遗传性疾病, 其病情的严重程度与ATP7B基因突变类型尚无明确定论，主要取决于铜在体内的蓄积量和积蓄时间，所以早期的诊疗对患者尤为重要^[5]。目前，对于转录调控如何影响肝豆状核变性的机制研究仍属于初期，本研究利用生物信息学分析方法筛选出7个在肝豆状核变性转录调控网络中位于核心地位的差异表达转录因子，并初步探讨了它们调控肝豆状核变性神经表型和肝脏表型可能的发生机制。目前，国内外关于肝豆状核变性的研究多集中在病例报道和突变基因型的确定，对于该病的转录调控机制和相关基因的功能研究还远远不够。未来将基于本研究结果，收集并纳入更多的肝豆状核变性就诊患者的血液和组织样本，进行大规模的高通量测序，结合生物信息学分析和分子生物学功能实验，多角度综合、详细地讨论肝豆状核变性患者的致病机制和转录调控机制。对于本研究发现的核心转录因子在高铜状态下具体如何对WD患者产生的影响，仍需要进一步通过实验进行系统地验证。