偏头痛皮质扩散性抑制的差异表达基因生物信息分析和治疗药物预测

周艳杰 张莉莉 杨 柳 王 月 肖哲曼

武汉大学人民医院，湖北 武汉 430060

偏头痛是神经系统常见的疾病，被世界卫生组织列为全球第六大致残性疾病［1］。最新的全球疾病负担调查（GBD 2016）显示，偏头痛的伤残损失健康生命年在神经系统疾病中位列第二［1］。皮质扩散性抑制（cortical spreading depression，CSD）是一种缓慢传播（2～5 mm/min）的神经元和胶质细胞的去极化波，以及随后引起的一段时间的皮层突触活动抑制。CSD一直被认为是偏头痛的主要病理生理机制之一，通常与偏头痛的先兆和偏头痛发作期有关。偏头痛CSD 发作引起的恶心、呕吐、畏光、畏声等先兆症状严重影响患者的工作、生活和学习［2］，因此寻找偏头痛CSD后候选生物标志物对偏头痛的早期诊断及发作期治疗极为重要。本研究通过提取美国国立生物中心的GEO 数据库中的相关数据，通过Deseq2 包以及clusterProfiler 包对家族性偏瘫性偏头痛1 型（familial hemiplegic migraine type 1，FHM1）小鼠诱发CSD 后皮质转录本进行分析，探讨偏头痛先兆及发作期的潜在机制［3-4］。此外，利用STRING 平台和CytoHubba 插件鉴定出10 个核心基因，并利用qPCR 对核心基因进行了进一步验证。通过对变化最显著的150个上调差异基因和43个下调差异基因进行CMap 分析，鉴定出了20 种可以逆转偏头痛差异基因的小分子化合物，为偏头痛的早期诊断及治疗提供了新的潜在药物靶点［5-6］。

1 材料和方法

1.1 数据来源数据来源于GEO 数据（https：//www.ncbi.nlm.nih.gov/geo），以“migraine”作为搜索词，获取到一个偏头痛小鼠模型数据集GSE67933［7］。该数据集 基 于Illumina GPL13112 平 台，由EISING 等［7］提交。本研究选择了其中6 个FHM1 偏头痛小鼠模型CSD 发作后的皮质样本和6 个野生型C57BL/6J 假手术后24 h的皮质样本，从样本中分离得到总RNA，再利用Illumina测序平台进行基因表达谱分析。

1.2 预处理数据与筛选差异表达基因基于GPL3112平台构建的信息，将Ensemble gene ID名称识别号转化为正式基因名称。利用R 软件中的sva包对样本进行背景表达值校正和数据归一化后，剔除掉离群样本，剩余3个FHM1偏头痛小鼠模型CSD发作后的皮质样本和3个野生型C57BL/6J假手术后24 h的皮质样本。再通过Deseq2包分析基因表达矩阵，采用Benjamini ＆ Hochberg错误发现率的方法对P 值进行校正，以降低假阳性率［8］。根据｜logFC｜＞1，Padjust＜0.05 的条件筛选出319 个差异表达基因（differential expression genes，DEGs），其中包括155个上调基因和35个下调基因。差异基因的可视化通过R软件中的Pheatmap包以及ggplot2包实现。

1.3 差异表达基因本体论分析和基因集富集分析基因本体论分析（gene ontology，GO）作为一种常用的生物信息学分析方法，可对基因产物及其功能特征进行注释，在分子生物学领域得到广泛应用。GO分析由生物过程（BP）、细胞成分（CC）、分子功能（MF）组成［9］。基因集富集分析（gene set enrichment analysis，GSEA）是一种用来计算预定的基因集在2组之间是否有统计学上显著的、一致性差异的方法，被广泛用于高通量测序样本的生物学功能分析［10］。本研究中的GO 和GSEA 分析均由R 软件中的clusterProfiler包实现［4］。

1.4 蛋白质-蛋白质相互作用网络分析DEGs中的蛋白质-蛋白质相互作用（protein-protein interaction，PPI）网络分析由STRING（The Search Tool for the Retrieval of Interacting Genes，STRING，version 11.0）数据库构建［11］。该数据库覆盖了约2 460 万个蛋白质，可以用来显示蛋白质之间的相互作用，从而探究蛋白质的直接和间接关联。以置信度得分≥0.4 和最小相互作用数≥1为筛选标准，得到PPI网络。核心基因使用Cytoscape 软件（Cytoscape 3.7.1）的插件Cytohubba进行识别，使用MCC算法，从PPI网络中挑选出与周围基因具有高度连通性前10个基因作为核心基因［5，12］。

1.5 动物及偏头痛模型的建立成年C57BL/6雄性小鼠（18～20 g）购于武汉大学人民医院实验动物中心。在无特定病原体（specific pathogen free，SPF）实验室条件下单独饲养，12 h/12 h 光照周期和正常饮食。适应1 周后，将小鼠随机分为对照组（n=3）和偏头痛组（n=3）。采用间歇性腹腔注射硝酸甘油（nitroglycerin，NTG）建立偏头痛小鼠模型，即每隔1天给予偏头痛组小鼠腹腔注射1次NTG 10 mg/kg，空白对照组腹腔注射等量的生理盐水，均注射5次。

1.6 小鼠皮层mRNA 表达检测用Trizol 试剂（Invitrogen，USA）从小鼠皮层组织中制备总RNA，并使 用First Strand cDNA Synthesis kit（Servicebio，China）制备cDNA。根据说明使用SYBR 预混EX Taq™Ⅱ试剂盒（Servicebio，China）进行实时定量PCR，并在Bio-Rad CFX Manager 2.1 实时PCR 系统（Bio-Rad，USA）中进行。使用采用2-ΔΔCt法计算相对基因表达量。引物序列见表1。

表1 rt-PCR引物序列Table 1 rt-PCR primer sequences

1.7 CMap分析Connectivity Map（CMap）（https：//portals.broadinstitute.org/cmap）是一个含有药物特异性基因表达谱的参考数据库，并将其与疾病特异性基因标记进行比较，可用于识别药物、基因与疾病之间的联系，从而确定潜在的治疗候选药物［6］。本研究中以P＜0.05 为筛选标准，Enrichment 评分按负值由高到低排序筛选出前3个化合物作为逆转DEGs表达改变的潜在靶点。

1.8 统计学分析DEGs的筛选以及GO和GSEA通路富集分析均在R 4.0.3 软件中执行。CMap 分析在CMap 网站中进行。rt-PCR 数据分析在GraphpadPrism 8.0 进行。计量资料以均数±标准误（±se）描述，以P＜0.05为差异有统计学意义。

2 结果

2.1 DEGs筛选结果与野生型假手术组小鼠相比，FHM1 型小鼠诱发CSD 24 h 后皮质样本中共筛选出190 个差异基因，其中上调差异基因有155 个，下调差异基因有35 个。从筛选出的所有DEGs中挑选出变化最显著的40 个差异基因绘制热图（图1），红色代表上调基因，蓝色代表下调基因；颜色越深，变化越显著。

图1 TOP 40 DEGs基因表达水平Figure 1 Gene expression levels of TOP 40 DEGs

2.2 DEGs的GO 分析 将筛选到的190 个差异基因进行GO 分析，如图2 所示，差异基因细胞组分显示，主要定位于细胞质、高尔基体、肌动蛋白细胞骨架以及神经元突触连接等。生物学过程显示，主要集中在细胞通讯、小分子GTP酶介导的信号转导、轴突的发育以及细胞内转运的调节等方面。分子功能显示，主要与转录共调节活动、蛋白质丝氨酸/苏氨酸激酶活性、酶激活活动以及小GTP 酶、Ras GTP 酶、ATP酶、泛素样蛋白转移酶等多种酶的活性有关。

图2 FHM1型小鼠诱发CSD后DEGs的GO分析Figure 2 GO analysis of DEGs after CSD induction in FHM1 mice

2.3 所有检测基因的KEGG富集分析运用GSEA分析鉴定FHM1 型小鼠诱发CSD 后和健康对照组小鼠皮质中差异有统计学意义的基因组（图3），结果显示，FHM1 型小鼠诱发CSD 后，基因主要富集在人类巨细胞病毒感染、EB 病毒感染、卡波西肉瘤病毒感染、朊病毒病等会引起机体免疫炎症反应的通路以及microRNA、细胞因子信号通路上，这些通路在FHM1型小鼠CSD诱发后被激活，现有待测基因并未富集到被抑制到的通路。

图3 GSEA富集分析水平Figure 3 GSEA enrichment analysis level

2.4 蛋白相互作用网络分析将差异基因分析得到的190个差异基因导入string平台并分析构建蛋白相互作用网络，得到188 个节点，199 条关系。将上述蛋白相互作用网络导入cytoscape 软件中，如图4 所示，红色代表上调，绿色代表下调。点的直径越大，P值越小。节点之间的连接颜色越红，代表两个节点之间的相互作用越强；反之，节点之间的连接颜色越绿，代表两个节点之间的相互作用越弱。

图4 DEGs蛋白相互作用网络Figure 4 DEGs-protein interaction network

2.5 核心基因识别CytoHubba中的MCC算法是已被证实的预测重要靶点较为精确的方法。为进一步鉴定DEGs 中的核心基因，利用Cytoscape 软件的CytoHubba 插件，按照MCC 方法计算并选取前10 个基因，并把这些基因作为核心基因，颜色越红，代表评分越高。

2.6 核心基因的鉴定为验证筛选出的核心基因是否在偏头痛中增加，构建偏头痛小鼠模型，并提取小鼠皮层的RNA。箱线图显示，相对于对照组小鼠模型，偏头痛小鼠模型皮质中9个中枢基因均显著变化，另外Mrps10 也有所增加，但差异无统计学意义（P＞0.05）。见图5。

图5 rt-PCR鉴定核心基因Figure 5 Identification of core genes by rt-PCR

2.7 CMap 分析为寻找潜在的小分子化合物来逆转DEGs 表达的改变，对上调最显著的前150个DEGs和所有下调的35个DEGs进行了CMap分析，筛选出对逆转偏头痛CSD发作最重要的20 种小分子化合物，其中排名靠前的3种小分子化合物：芹菜素、1，4-屈烯醌、非诺特罗。

3 讨论

本研究利用已报道的GEO 芯片数据集GSE67933 进行分析，并通过提取其中的6只FHM1型小鼠诱发CSD后和6只健康对照小鼠皮质的高通量测序样本进行差异基因的分析。去批次和背景归一化后删除掉离群样本，剩余样本通过R软件的Deseq2包筛选出了190个差异基因，其中包括155 个上调基因和35 个下调基因。对差异基因进行GO 分析显示，主要定位于细胞质、高尔基体、肌动蛋白细胞骨架以及神经元突触连接等。

GO 分析的生物学过程显示，主要集中在细胞通讯、小分子GTP酶介导的信号转导、轴突的发育以及细胞内转运的调节等方面。分子功能显示，主要与转录共调节活动、蛋白质丝氨酸/苏氨酸激酶活性、酶激活活动以及小GTP 酶、Ras GTP 酶、ATP酶、泛素样蛋白转移酶等多种酶的活性有关。GSEA分析表明FHM1型小鼠诱发CSD后皮质DEGs主要富集在一些会引起机体免疫炎症反应的通路以及microRNA、细胞因子信号通路上。这些通路在FHM1型小鼠CSD诱发后被激活，现有待测基因并未富集到FHM1 型小鼠CSD发作后被抑制的通路。

通过构建DEGs 的蛋白相互作用网络，并计算MCC 评分得出10 个与其他节点作用最多的基因并将这些基因作为核心 基 因，分 别 是Rplp0、Rpl4、Eef1b2、Eif5a、Mrps10、Rpl17、Rpl15、Tsfm、Src、Hsp90ab1。通过 人 类 基 因 数 据 库GeneCards （https：//www.genecards.org/）对这10 个核心基因注释发现这些基因主要与RNA 的结合、核糖体途径以及蛋白质的翻译过程有关［13］。其中Hsp90ab1 即热休克蛋白90 Alpha 家族B 类成员1，通过编码细胞溶质90 kDa 热休克蛋白（Hsp90）在多种疾病中发挥作用，在疼痛领域中已得到广泛研究。最近的一项实验表明抑制脊髓中的Hsp90 后ERK-RSK 通路被激活，进而增强了吗啡的镇痛作用［14］。一项关于单关节炎大鼠模型的实验表明，Hsp90 抑制剂可在最初几个小时内减轻MA引起的机械性异常性疼痛［15］，表明Hsp90ab1有可能与偏头痛的痛觉敏化有关。Src编码酪氨酸蛋白激酶，与神经病理性疼痛、慢性疼痛、癌痛多种疼痛反应有关［16-18］。许多研究表明，酪氨酸蛋白激酶的激活或磷酸化在疼痛超敏反应中起重要作用［19］。Rplp0、Rp14、Rpl17和Rpl15主要参与核糖体结构和蛋白质合成。Eefib2编码一种翻译延伸因子，称为真核翻译延伸因子1beta2，参与将氨酰化tRNA 转移到核糖体。Mrps10编码哺乳动物线粒体核糖体蛋白的核基因，有助于线粒体内的蛋白质合成，与C 反应蛋白、中性粒细胞、白细胞等多种炎症反应有关。本研究用偏头痛小鼠模型对Mrps10验证未得到有统计学意义的增加，可能是样本量较小的原因。Eif5a 参与翻译延伸的mRNA 结合蛋白，在mRNA 更新水平上具有重要功能，参与肌动蛋白动力学和细胞周期进程，并可能参与应激反应和维持细胞壁完整性的途径，其与合成蛋白Sdcbp 一起，作为p53/TP53 和p53/TP53 依赖性细胞凋亡的调节剂起作用，还调节TNF-α介导的细胞凋亡，介导多胺对神经元过程延伸和存活的影响，可能在大脑发育和功能以及骨骼肌干细胞分化中发挥重要作用。Tsfm编码线粒体翻译延伸因子，在线粒体蛋白翻译的延伸步骤中，催化翻译延伸因子Tu 上鸟嘌呤核苷酸的交换，研究表明Tsfm基因突变会引起多动症、脑心肌病、儿童共济失调以及舞蹈病等多种神经系统疾病［20］。结合GO 和GSEA分析可以推断偏头痛CSD发作后关键基因的上调引起下游免疫细胞的蛋白质合成，激活了炎症反应，从而进一步在偏头痛疼痛反应的启动、维持和延伸中发挥作用。

为进一步探究与偏头痛CSD发作有关的偏头痛先兆以及急性偏头痛的潜在治疗方案，进行CMap分析，通过对变化最显著的前5个上调DEGs和35个下调DEGs 进行预测，得到最有可能逆转偏头痛CSD 3 种小分子化合物是芹菜素、1,4-屈烯醌、非诺特罗。

芹菜素是一种黄酮类物质，具有抑制血管平滑肌细胞的增殖、抗氧化、抑制肿瘤细胞的生长和转移、镇静、改善睡眠等作用［21］，目前主要运用于预防高血压、动脉硬化、心脑血管疾病等方面。最近许多研究表明芹菜素对偏头痛有一定的镇痛作用［22］。一项关于偏头痛的植物疗法研究显示，使用以芹菜素为主要成分的洋甘菊油凝胶30 min 后，患者的疼痛、恶心、呕吐、畏光和畏声等现象显著降低，表明芹菜素对偏头痛有一定的治疗功效［22］。另一项研究表明，芹菜素主要通过以下机制缓解偏头痛：（1）抑制活化巨噬细胞中iNOS 的表达，导致NO 释放和合成受阻。NO 在刺激中枢敏化方面具有重要作用，其合成的阻断或减少可以帮助治疗偏头痛发作和疼痛。（2）选择性抑制内源性前列腺素E2（PGE2）水平，阻断中枢神经元和外周脑膜伤害性感受器的致敏作用。（3）具有与皮质类固醇（如氢化可的松）一样强的抗炎作用，通过缓解作用部位的炎症减轻偏头痛［23］。

1,4-屈烯醌（1,4-CQ）是芳烃受体（aryl hydrocarbon receptor，AhR）的激活剂。AhR是一种配体激活的转录因子，具有许多生理功能，包括调节胃肠道和调节肠脑轴。许多研究表明AhR 的活化与偏头痛的发病机制有关［24］。最近一项关于偏头痛和功能性胃肠疾病中色氨酸代谢的研究表明，偏头痛中存在犬尿氨酸途径代谢的变化，而犬尿氨酸途径的产物是AhR 的重要配体，犬尿氨酸类似物可竞争性拮抗NMDA受体，减少偏头痛中神经元的过度活跃，从而消除硝酸甘油诱导的痛觉过敏［25］。此外，AhR 还可通过介导毒素反应调节T细胞免疫，而偏头痛的发病机制一直被认为与免疫炎症有关［26］，表明1,4-屈烯醌有望成为偏头痛治疗的候选药物靶点之一。

非诺特罗是一种选择性β2受体激动剂，属间羟异丙肾上腺素的衍生物，临床上主要用于支气管哮喘的治疗。目前有研究发现非诺特罗等β2-肾上腺素能受体（β2-AR）能在很大程度上缓解神经病理性疼痛小鼠模型的疼痛反应［27］，认为β2-肾上腺素能受体刺激物可作为潜在的神经性疼痛的新疗法。由此可以推测非诺特罗或许在偏头痛中也有一定的治疗效果。

本研究通过生物信息学分析确定了10个FHM1型小鼠CSD发作后的10个核心基因以及偏头痛CSD发作的潜在机制，并预测了可用来逆转DEGs 的20种小分子化合物，但仍需要进一步的研究来验证核心基因和深入的致病机制。本研究可为与CSD发作有关的先兆偏头痛及急性偏头痛的机制和治疗提供潜在的靶点。