偏頭痛皮質擴散性抑制的差異表達基因生物信息分析和治療藥物預測

周艷杰 張莉莉 楊 柳 王 月 肖哲曼

武漢大學人民醫(yī)院，湖北 武漢 430060

偏頭痛是神經(jīng)系統(tǒng)常見的疾病，被世界衛(wèi)生組織列為全球第六大致殘性疾病［1］。最新的全球疾病負擔調查（GBD 2016）顯示，偏頭痛的傷殘損失健康生命年在神經(jīng)系統(tǒng)疾病中位列第二［1］。皮質擴散性抑制（cortical spreading depression，CSD）是一種緩慢傳播（2～5 mm/min）的神經(jīng)元和膠質細胞的去極化波，以及隨后引起的一段時間的皮層突觸活動抑制。CSD一直被認為是偏頭痛的主要病理生理機制之一，通常與偏頭痛的先兆和偏頭痛發(fā)作期有關。偏頭痛CSD 發(fā)作引起的惡心、嘔吐、畏光、畏聲等先兆癥狀嚴重影響患者的工作、生活和學習［2］，因此尋找偏頭痛CSD后候選生物標志物對偏頭痛的早期診斷及發(fā)作期治療極為重要。本研究通過提取美國國立生物中心的GEO 數(shù)據(jù)庫中的相關數(shù)據(jù)，通過Deseq2 包以及clusterProfiler 包對家族性偏癱性偏頭痛1 型（familial hemiplegic migraine type 1，F(xiàn)HM1）小鼠誘發(fā)CSD 后皮質轉錄本進行分析，探討偏頭痛先兆及發(fā)作期的潛在機制［3-4］。此外，利用STRING 平臺和CytoHubba 插件鑒定出10 個核心基因，并利用qPCR 對核心基因進行了進一步驗證。通過對變化最顯著的150個上調差異基因和43個下調差異基因進行CMap 分析，鑒定出了20 種可以逆轉偏頭痛差異基因的小分子化合物，為偏頭痛的早期診斷及治療提供了新的潛在藥物靶點［5-6］。

1 材料和方法

1.1 數(shù)據(jù)來源數(shù)據(jù)來源于GEO 數(shù)據(jù)（https：//www.ncbi.nlm.nih.gov/geo），以“migraine”作為搜索詞，獲取到一個偏頭痛小鼠模型數(shù)據(jù)集GSE67933［7］。該數(shù)據(jù)集 基 于Illumina GPL13112 平 臺，由EISING 等［7］提交。本研究選擇了其中6 個FHM1 偏頭痛小鼠模型CSD 發(fā)作后的皮質樣本和6 個野生型C57BL/6J 假手術后24 h的皮質樣本，從樣本中分離得到總RNA，再利用Illumina測序平臺進行基因表達譜分析。

1.2 預處理數(shù)據(jù)與篩選差異表達基因基于GPL3112平臺構建的信息，將Ensemble gene ID名稱識別號轉化為正式基因名稱。利用R 軟件中的sva包對樣本進行背景表達值校正和數(shù)據(jù)歸一化后，剔除掉離群樣本，剩余3個FHM1偏頭痛小鼠模型CSD發(fā)作后的皮質樣本和3個野生型C57BL/6J假手術后24 h的皮質樣本。再通過Deseq2包分析基因表達矩陣，采用Benjamini ＆ Hochberg錯誤發(fā)現(xiàn)率的方法對P 值進行校正，以降低假陽性率［8］。根據(jù)｜logFC｜＞1，Padjust＜0.05 的條件篩選出319 個差異表達基因（differential expression genes，DEGs），其中包括155個上調基因和35個下調基因。差異基因的可視化通過R軟件中的Pheatmap包以及ggplot2包實現(xiàn)。

1.3 差異表達基因本體論分析和基因集富集分析基因本體論分析（gene ontology，GO）作為一種常用的生物信息學分析方法，可對基因產(chǎn)物及其功能特征進行注釋，在分子生物學領域得到廣泛應用。GO分析由生物過程（BP）、細胞成分（CC）、分子功能（MF）組成［9］。基因集富集分析（gene set enrichment analysis，GSEA）是一種用來計算預定的基因集在2組之間是否有統(tǒng)計學上顯著的、一致性差異的方法，被廣泛用于高通量測序樣本的生物學功能分析［10］。本研究中的GO 和GSEA 分析均由R 軟件中的clusterProfiler包實現(xiàn)［4］。

1.4 蛋白質-蛋白質相互作用網(wǎng)絡分析DEGs中的蛋白質-蛋白質相互作用（protein-protein interaction，PPI）網(wǎng)絡分析由STRING（The Search Tool for the Retrieval of Interacting Genes，STRING，version 11.0）數(shù)據(jù)庫構建［11］。該數(shù)據(jù)庫覆蓋了約2 460 萬個蛋白質，可以用來顯示蛋白質之間的相互作用，從而探究蛋白質的直接和間接關聯(lián)。以置信度得分≥0.4 和最小相互作用數(shù)≥1為篩選標準，得到PPI網(wǎng)絡。核心基因使用Cytoscape 軟件（Cytoscape 3.7.1）的插件Cytohubba進行識別，使用MCC算法，從PPI網(wǎng)絡中挑選出與周圍基因具有高度連通性前10個基因作為核心基因［5，12］。

1.5 動物及偏頭痛模型的建立成年C57BL/6雄性小鼠（18～20 g）購于武漢大學人民醫(yī)院實驗動物中心。在無特定病原體（specific pathogen free，SPF）實驗室條件下單獨飼養(yǎng)，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。根據(jù)說明使用SYBR 預混EX Taq?Ⅱ試劑盒（Servicebio，China）進行實時定量PCR，并在Bio-Rad CFX Manager 2.1 實時PCR 系統(tǒng)（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）是一個含有藥物特異性基因表達譜的參考數(shù)據(jù)庫，并將其與疾病特異性基因標記進行比較，可用于識別藥物、基因與疾病之間的聯(lián)系，從而確定潛在的治療候選藥物［6］。本研究中以P＜0.05 為篩選標準，Enrichment 評分按負值由高到低排序篩選出前3個化合物作為逆轉DEGs表達改變的潛在靶點。

1.8 統(tǒng)計學分析DEGs的篩選以及GO和GSEA通路富集分析均在R 4.0.3 軟件中執(zhí)行。CMap 分析在CMap 網(wǎng)站中進行。rt-PCR 數(shù)據(jù)分析在GraphpadPrism 8.0 進行。計量資料以均數(shù)±標準誤（±se）描述，以P＜0.05為差異有統(tǒng)計學意義。

2 結果

2.1 DEGs篩選結果與野生型假手術組小鼠相比，F(xiàn)HM1 型小鼠誘發(fā)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 所示，差異基因細胞組分顯示，主要定位于細胞質、高爾基體、肌動蛋白細胞骨架以及神經(jīng)元突觸連接等。生物學過程顯示，主要集中在細胞通訊、小分子GTP酶介導的信號轉導、軸突的發(fā)育以及細胞內轉運的調節(jié)等方面。分子功能顯示，主要與轉錄共調節(jié)活動、蛋白質絲氨酸/蘇氨酸激酶活性、酶激活活動以及小GTP 酶、Ras GTP 酶、ATP酶、泛素樣蛋白轉移酶等多種酶的活性有關。

圖2 FHM1型小鼠誘發(fā)CSD后DEGs的GO分析Figure 2 GO analysis of DEGs after CSD induction in FHM1 mice

2.3 所有檢測基因的KEGG富集分析運用GSEA分析鑒定FHM1 型小鼠誘發(fā)CSD 后和健康對照組小鼠皮質中差異有統(tǒng)計學意義的基因組（圖3），結果顯示，F(xiàn)HM1 型小鼠誘發(fā)CSD 后，基因主要富集在人類巨細胞病毒感染、EB 病毒感染、卡波西肉瘤病毒感染、朊病毒病等會引起機體免疫炎癥反應的通路以及microRNA、細胞因子信號通路上，這些通路在FHM1型小鼠CSD誘發(fā)后被激活，現(xiàn)有待測基因并未富集到被抑制到的通路。

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

2.4 蛋白相互作用網(wǎng)絡分析將差異基因分析得到的190個差異基因導入string平臺并分析構建蛋白相互作用網(wǎng)絡，得到188 個節(jié)點，199 條關系。將上述蛋白相互作用網(wǎng)絡導入cytoscape 軟件中，如圖4 所示，紅色代表上調，綠色代表下調。點的直徑越大，P值越小。節(jié)點之間的連接顏色越紅，代表兩個節(jié)點之間的相互作用越強；反之，節(jié)點之間的連接顏色越綠，代表兩個節(jié)點之間的相互作用越弱。

圖4 DEGs蛋白相互作用網(wǎng)絡Figure 4 DEGs-protein interaction network

2.5 核心基因識別CytoHubba中的MCC算法是已被證實的預測重要靶點較為精確的方法。為進一步鑒定DEGs 中的核心基因，利用Cytoscape 軟件的CytoHubba 插件，按照MCC 方法計算并選取前10 個基因，并把這些基因作為核心基因，顏色越紅，代表評分越高。

2.6 核心基因的鑒定為驗證篩選出的核心基因是否在偏頭痛中增加，構建偏頭痛小鼠模型，并提取小鼠皮層的RNA。箱線圖顯示，相對于對照組小鼠模型，偏頭痛小鼠模型皮質中9個中樞基因均顯著變化，另外Mrps10 也有所增加，但差異無統(tǒng)計學意義（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發(fā)作最重要的20 種小分子化合物，其中排名靠前的3種小分子化合物：芹菜素、1，4-屈烯醌、非諾特羅。

3 討論

本研究利用已報道的GEO 芯片數(shù)據(jù)集GSE67933 進行分析，并通過提取其中的6只FHM1型小鼠誘發(fā)CSD后和6只健康對照小鼠皮質的高通量測序樣本進行差異基因的分析。去批次和背景歸一化后刪除掉離群樣本，剩余樣本通過R軟件的Deseq2包篩選出了190個差異基因，其中包括155 個上調基因和35 個下調基因。對差異基因進行GO 分析顯示，主要定位于細胞質、高爾基體、肌動蛋白細胞骨架以及神經(jīng)元突觸連接等。

GO 分析的生物學過程顯示，主要集中在細胞通訊、小分子GTP酶介導的信號轉導、軸突的發(fā)育以及細胞內轉運的調節(jié)等方面。分子功能顯示，主要與轉錄共調節(jié)活動、蛋白質絲氨酸/蘇氨酸激酶活性、酶激活活動以及小GTP 酶、Ras GTP 酶、ATP酶、泛素樣蛋白轉移酶等多種酶的活性有關。GSEA分析表明FHM1型小鼠誘發(fā)CSD后皮質DEGs主要富集在一些會引起機體免疫炎癥反應的通路以及microRNA、細胞因子信號通路上。這些通路在FHM1型小鼠CSD誘發(fā)后被激活，現(xiàn)有待測基因并未富集到FHM1 型小鼠CSD發(fā)作后被抑制的通路。

通過構建DEGs 的蛋白相互作用網(wǎng)絡，并計算MCC 評分得出10 個與其他節(jié)點作用最多的基因并將這些基因作為核心 基 因，分 別 是Rplp0、Rpl4、Eef1b2、Eif5a、Mrps10、Rpl17、Rpl15、Tsfm、Src、Hsp90ab1。通過 人 類 基 因 數(shù) 據(jù) 庫GeneCards （https：//www.genecards.org/）對這10 個核心基因注釋發(fā)現(xiàn)這些基因主要與RNA 的結合、核糖體途徑以及蛋白質的翻譯過程有關［13］。其中Hsp90ab1 即熱休克蛋白90 Alpha 家族B 類成員1，通過編碼細胞溶質90 kDa 熱休克蛋白（Hsp90）在多種疾病中發(fā)揮作用，在疼痛領域中已得到廣泛研究。最近的一項實驗表明抑制脊髓中的Hsp90 后ERK-RSK 通路被激活，進而增強了嗎啡的鎮(zhèn)痛作用［14］。一項關于單關節(jié)炎大鼠模型的實驗表明，Hsp90 抑制劑可在最初幾個小時內減輕MA引起的機械性異常性疼痛［15］，表明Hsp90ab1有可能與偏頭痛的痛覺敏化有關。Src編碼酪氨酸蛋白激酶，與神經(jīng)病理性疼痛、慢性疼痛、癌痛多種疼痛反應有關［16-18］。許多研究表明，酪氨酸蛋白激酶的激活或磷酸化在疼痛超敏反應中起重要作用［19］。Rplp0、Rp14、Rpl17和Rpl15主要參與核糖體結構和蛋白質合成。Eefib2編碼一種翻譯延伸因子，稱為真核翻譯延伸因子1beta2，參與將氨?；痶RNA 轉移到核糖體。Mrps10編碼哺乳動物線粒體核糖體蛋白的核基因，有助于線粒體內的蛋白質合成，與C 反應蛋白、中性粒細胞、白細胞等多種炎癥反應有關。本研究用偏頭痛小鼠模型對Mrps10驗證未得到有統(tǒng)計學意義的增加，可能是樣本量較小的原因。Eif5a 參與翻譯延伸的mRNA 結合蛋白，在mRNA 更新水平上具有重要功能，參與肌動蛋白動力學和細胞周期進程，并可能參與應激反應和維持細胞壁完整性的途徑，其與合成蛋白Sdcbp 一起，作為p53/TP53 和p53/TP53 依賴性細胞凋亡的調節(jié)劑起作用，還調節(jié)TNF-α介導的細胞凋亡，介導多胺對神經(jīng)元過程延伸和存活的影響，可能在大腦發(fā)育和功能以及骨骼肌干細胞分化中發(fā)揮重要作用。Tsfm編碼線粒體翻譯延伸因子，在線粒體蛋白翻譯的延伸步驟中，催化翻譯延伸因子Tu 上鳥嘌呤核苷酸的交換，研究表明Tsfm基因突變會引起多動癥、腦心肌病、兒童共濟失調以及舞蹈病等多種神經(jīng)系統(tǒng)疾病［20］。結合GO 和GSEA分析可以推斷偏頭痛CSD發(fā)作后關鍵基因的上調引起下游免疫細胞的蛋白質合成，激活了炎癥反應，從而進一步在偏頭痛疼痛反應的啟動、維持和延伸中發(fā)揮作用。

為進一步探究與偏頭痛CSD發(fā)作有關的偏頭痛先兆以及急性偏頭痛的潛在治療方案，進行CMap分析，通過對變化最顯著的前5個上調DEGs和35個下調DEGs 進行預測，得到最有可能逆轉偏頭痛CSD 3 種小分子化合物是芹菜素、1,4-屈烯醌、非諾特羅。

芹菜素是一種黃酮類物質，具有抑制血管平滑肌細胞的增殖、抗氧化、抑制腫瘤細胞的生長和轉移、鎮(zhèn)靜、改善睡眠等作用［21］，目前主要運用于預防高血壓、動脈硬化、心腦血管疾病等方面。最近許多研究表明芹菜素對偏頭痛有一定的鎮(zhèn)痛作用［22］。一項關于偏頭痛的植物療法研究顯示，使用以芹菜素為主要成分的洋甘菊油凝膠30 min 后，患者的疼痛、惡心、嘔吐、畏光和畏聲等現(xiàn)象顯著降低，表明芹菜素對偏頭痛有一定的治療功效［22］。另一項研究表明，芹菜素主要通過以下機制緩解偏頭痛：（1）抑制活化巨噬細胞中iNOS 的表達，導致NO 釋放和合成受阻。NO 在刺激中樞敏化方面具有重要作用，其合成的阻斷或減少可以幫助治療偏頭痛發(fā)作和疼痛。（2）選擇性抑制內源性前列腺素E2（PGE2）水平，阻斷中樞神經(jīng)元和外周腦膜傷害性感受器的致敏作用。（3）具有與皮質類固醇（如氫化可的松）一樣強的抗炎作用，通過緩解作用部位的炎癥減輕偏頭痛［23］。

1,4-屈烯醌（1,4-CQ）是芳烴受體（aryl hydrocarbon receptor，AhR）的激活劑。AhR是一種配體激活的轉錄因子，具有許多生理功能，包括調節(jié)胃腸道和調節(jié)腸腦軸。許多研究表明AhR 的活化與偏頭痛的發(fā)病機制有關［24］。最近一項關于偏頭痛和功能性胃腸疾病中色氨酸代謝的研究表明，偏頭痛中存在犬尿氨酸途徑代謝的變化，而犬尿氨酸途徑的產(chǎn)物是AhR 的重要配體，犬尿氨酸類似物可競爭性拮抗NMDA受體，減少偏頭痛中神經(jīng)元的過度活躍，從而消除硝酸甘油誘導的痛覺過敏［25］。此外，AhR 還可通過介導毒素反應調節(jié)T細胞免疫，而偏頭痛的發(fā)病機制一直被認為與免疫炎癥有關［26］，表明1,4-屈烯醌有望成為偏頭痛治療的候選藥物靶點之一。

非諾特羅是一種選擇性β2受體激動劑，屬間羥異丙腎上腺素的衍生物，臨床上主要用于支氣管哮喘的治療。目前有研究發(fā)現(xiàn)非諾特羅等β2-腎上腺素能受體（β2-AR）能在很大程度上緩解神經(jīng)病理性疼痛小鼠模型的疼痛反應［27］，認為β2-腎上腺素能受體刺激物可作為潛在的神經(jīng)性疼痛的新療法。由此可以推測非諾特羅或許在偏頭痛中也有一定的治療效果。

本研究通過生物信息學分析確定了10個FHM1型小鼠CSD發(fā)作后的10個核心基因以及偏頭痛CSD發(fā)作的潛在機制，并預測了可用來逆轉DEGs 的20種小分子化合物，但仍需要進一步的研究來驗證核心基因和深入的致病機制。本研究可為與CSD發(fā)作有關的先兆偏頭痛及急性偏頭痛的機制和治療提供潛在的靶點。