慢性疼痛導致抑郁的腦區基因共表達網絡分析 *

錢召強 李永豐 劉一輝 任 維 田英芳△

（1 陜西師范大學生命科學學院，西安710119；2 陜西中醫藥大學針灸推拿學院，西安712046；3 陜西師范大學現代教學技術教育部重點實驗室，西安710062）

神經病理性疼痛 (neuropathic pain, NP) 是一種由軀體感覺神經系統的損傷或疾病引起的慢性疼痛，臨床表現主要為自發性疼痛、痛覺過敏和痛覺超敏[1]。神經病理性疼痛在人群中的患病率高達7%～10%，給個人和社會帶來了沉重的經濟負擔[2]。該疾病不僅導致病患生活質量下降、健康不佳等，還明顯導致了病患的情緒低落，約有50%的病人發生焦慮和抑郁[3]。同時，焦慮和抑郁等負性情緒加劇了痛覺感受。目前，慢性疼痛誘發焦慮和抑郁的神經機制尚不明確。有研究表明慢性疼痛與情緒、動機、獎賞相關的幾個腦區的適應性有關。伏隔核(nucleus accumbens, NAc) 位于基底核與邊緣系統交界處，是獎賞通路的重要組成部分，參與獎賞、自發活動等行為，同時在痛覺信息調制中具有重要作用[4,5]。內側前額葉皮層 (medial prefrontal cortex, mPFC)與工作記憶、決策、情緒調控和注意力等高級認知功能有關，并且參與了痛覺信息處理[5,6]。中腦導水管周圍灰質 (periaqueductal grey, PAG) 在傷害性感覺的下行調節中起著重要作用，并調節情緒行為[7]。

基因芯片和轉錄組高通量測序技術 (RNA Sequencing, RNA-seq) 可以大規模地檢測基因表達譜，因此已經廣泛地應用于各種疾病的分子機制研究。Descalzi 等[8]研究了慢性疼痛導致抑郁癥的NAc、PAG 及mPFC 等三個腦區的轉錄組及模式，為相關研究提供了契機。RNA-seq 涉及到大量的基因數據，傳統的基于兩兩比較獲得差異表達基因的方法不能有效處理海量數據，因而引入了網絡分析，其中加權基因共表達網絡分析 (weighted gene coexpression network analysis, WGCNA) 被廣泛應用。WGCNA是一種系統生物學研究方法，可識別不同樣本間表達變化趨勢相似的基因，是大規模數據集中進行基因及疾病預測分析的有效策略。該方法充分利用了轉錄基因的信息，探究基因功能及基因間的相互作用關系，因而更接近生物體內的真實情況，能發現某些表達差異不顯著但有重要生物學意義的基因，可以獲得更為有意義的結果。因此，本研究使用Descalzi 等公開的RNA-seq 數據進行WGCNA 分析，研究神經病理性疼痛導致抑郁的小鼠NAC、PAG及mPFC 三個腦區的基因表達特征，以期對慢性疼痛導致抑郁的機制及治療提供新思路。

方 法

1. 數據來源和動物模型

本研究使用的轉錄組數據下載自美國國家生物技術信息中心 (NCBI) 的基因表達數據庫 (gene expression omnibus, GEO) (www.ncbi.nlm.nih.gov/geo)，存儲號為GSE91396，由Descalzi 等分享[8]，為基于Illumina HiSeq 2000 平臺的高通量測序技術獲得的RNA-seq 數據，共包含36 個樣本的數據。

實驗采用8～9 周齡12 只雄性成年C57BL/6J小鼠構建坐骨神經分支損傷 (spared nerve injury,SNI) 模型（6只為假手術對照組，6只為SNI手術組）。手術時用三溴乙醇進行全身麻醉，在小鼠左后腿大腿中部切開皮膚和肌肉，暴露坐骨神經及其三個分支。用6.0 縫合線結扎腓總神經和腓腸總神經，剪斷并切除1～2 mm 神經纖維，脛神經保持完好無損，用4.0 縫合線縫合皮膚。假手術的小鼠接受同樣的操作程序，但所有的神經都完好無損。術后第14、28、42、63 天時檢測痛覺閾值，SNI 手術組的閾值維持在較低水平，與對照組相比有顯著性差異，術后第9 周進行強迫游泳、糖水偏愛等實驗檢測，SNI 手術組表現出抑郁樣行為，第77 天處死小鼠分離兩側NAc、PAG 及mPFC，使用TRIzol 法分離總RNA，建庫后進行RNA 測序，使用TopHat2 軟件將測序片段與小鼠基因組序列比對獲得基因表達數量，用HTSeq 軟件進行基因注釋。

2. 加權共表達基因網絡構建

使用R 軟件的WGCNA 包（版本1.69）進行加權共表達網絡分析[9]。首先使用hclust 函數進行樣本聚類，剔除離群樣本。隨后用pickSoftThreshold 函數計算構建網絡的軟閾值 (power)，無尺度網絡擬合指數R2≥0.8，軟閾值使鄰接矩陣為0～1 之間的連續值，從而使構建的網絡符合冪分布，更接近實際的生物網絡狀態。然后用blockwiseModules函數構建拓撲重疊矩陣 (topological overlap matrix,TOM)，使用動態樹剪切算法進行聚類，將表達模式相似的基因分組到同一共表達模塊，每種模塊分配不同的顏色進行區分。其中，模塊最少基因數為30，合并相似模塊的閾值為0.25，其余使用默認參數。共表達模塊中的基因具有很高的連通性，同一模塊中的基因可能具有相似的生物學功能。

3. 識別與腦區相關性高的模塊

NAc、PAG 及mPFC 三個腦區在慢性疼痛導致的抑郁癥中發揮不同作用，為了確定與三個腦區顯著相關的基因模塊，對每個模塊進行主成分分析(principle component analysis, PCA)，將模塊特征基因值 (module eigengene, ME) 即模塊的第一主成分[10]與表型性狀（即3 個腦區）相關聯，與性狀相關性最高（R＞ 0.95 且P≤0.05）的模塊即為組織特異性模塊。

4. 模塊基因功能和通路富集分析

將與腦區顯著相關的基因模塊進行模塊功能和通路富集分析，以進一步探索各模塊的生物學功能。本研究使用R 軟件的clusterProfiler 包對各感興趣模塊進行GO (gene ontology) 和KEGG (kyoto encyclopedia of genes and genomes) 分析，以P≤0.05 為顯著性閾值。GO 分析在生物過程(biological process,BP)、分子功能 (molecular function, MF)、細胞成分(cellular component, CC) 等3 個方面進行基因及其產物功能的注釋；KEGG 數據庫集合了基因、蛋白質、化學組分等信息及其相互作用、反應和關系網絡，對基因功能和代謝通路進行注釋分析。在結果中分別選擇排名前10 個條目 (term) 進行分析。

5. 篩選樞紐基因

在模塊內拓撲指數最高也即與其他基因的連接數量最大基因為樞紐基因，在模塊內的連通性最高，因而有比較重要的作用。使用Cytoscape 中cytoHubba 插件的MCC 算法，分析WGCNA 得到的網絡數據，取連通性前10 的基因為樞紐基因。

結 果

1. 加權共表達基因網絡構建

對下載的數據進行處理，刪除基因表達量為0 的數據，共獲得23, 503 個基因。36 個數據經計算數據相關系數后聚類，刪除聚類異常的1 個數據（編號為sample19），最終用于分析的有35 個數據。根據絕對中位差選取前8000 個表達量最大的基因，進行下一步WGCNA 分析。根據無尺度網絡的擬合度0.8，得到構建網絡的軟閾值為16，平均連通性趨近于0（見圖1）。通過動態剪接法獲得9 個模塊（見圖2），對應的顏色分別為black、blue、brown、green、grey、pink、red、turquoise、yellow，模塊中的基因數在96～3813 不等，其中灰色模塊為不能聚類到其他模塊的基因，因而后續分析中不進行分析。

圖1 網絡軟閾值和連通性Fig. 1 Soft-thresholdingpowers and mean connectivity

圖2 基因聚類樹狀圖以及模塊Fig. 2 Clustering dendrogram of genes, with assigned module colors

2. 識別與腦區相關性高的模塊

為了確定每個模塊的重要性，將所獲得的模塊與三個腦區進行相關性分析（見圖3）。從熱圖中可以看出brown模塊（1029個基因）與NAc (R= 0.96,P= 3e-19)、turquoise 模塊（3813 個基因）與PAG(R= 0.97,P= 4e-22)、yellow 模塊（834 個基因）與mPFC (R= 0.99,P= 6e-27) 高度正相關，提示三個腦區有著不同的基因表達模式，以此三個模塊為目標模塊做進一步分析。

圖3 基因共表達網絡模塊與腦區的相關性熱圖括號內數字表示P 值，括號外數字表示相關系數R 值Fig. 3 Heatmap of the correlation between modules and 3 brain regions Number in and outside the bracket represents P value and Pearson coefficient R value, respectively.

3. 模塊功能和通路富集分析

為了得到三個待研究模塊的生物學功能，用R軟件的ClusterProfile 包進行GO 和KEGG 分析。GO 富集分析結果表明（見圖4），NAc 腦區在認知、突觸膜電位、神經遞質轉運蛋白、學習記憶、囊泡介導的轉運蛋白調節等生物學過程中發揮重要作用，具有陽離子通道活性、底物特異性通道活性、門控通道活性、小GTP 酶結合、Ras GTP 酶結合等分子功能，主要的細胞組分為不對稱突觸、突觸后膜、突觸后特化、突觸后密度等；PAG 腦區在突觸組織、細胞投射組織的正向調節、樹突發育、有機化合物氧化產生的能量、磷酸核糖代謝過程等生物學過程中發揮重要作用，具有泛素蛋白連接酶結合、小GTP 酶結合、突觸的結構成分、Ras GTP 酶結合、突觸后的結構成分、磷脂結合等分子功能，主要細胞組分為髓鞘、不對稱突觸、突觸后致密區、突觸后分化、細胞器內膜、突觸后膜等；mPFC 腦區在軸突生成、突觸組織、神經遞質轉運、突觸小泡循環的調節、小泡介導的突觸轉運等生物學過程中發揮重要作用，具有GTP 酶激活因子活性、小GTP酶結合、Ras GTP酶結合、離子通道結合等分子功能，主要細胞組分有胞外囊泡、突觸囊泡、轉運囊泡膜、突觸囊泡膜、突觸前囊泡膜等。三個腦區都分析得到小GTP 酶結合 (GO: 0031267)、Ras GTP 酶結合(GO: 0017016) 兩個分子功能，提示在此模型中可能有重要作用。

圖4 NAc、mPFC 和PAG 腦區的GO 功能注釋Fig. 4 GO functional analysis for genes in NAc, mPFC and PAG, respectively

KEGG 富集結果表明（見圖5），NAc 腦區的基因顯著富集在長時程增強、嗎啡成癮、膽堿能突觸、谷氨酸能突觸、多巴胺能突觸等通路；PAG 腦區的基因顯著富集在帕金森病、亨廷頓病、內質網蛋白加工等通路；mPFC 腦區的基因顯著富集在逆向內源性大麻素信號、突觸囊泡周期、谷氨酸能突觸、GABA 能突觸、多巴胺能突觸等通路。

圖5 NAc、mPFC 和PAG 腦區的KEGG 富集分析Fig. 5 KEGG pathway enrich analysis for genes in NAc,mPFC and PAG, respectively

4. 篩選樞紐基因

分析共表達基因相互作用網絡中連接數，用Cytoscape 軟件篩選連通性最高的前10 個基因為樞紐基因。NAc 腦區的樞紐基因為Drd1a、Strip2、Slmap、Gpr88、Cacna2d3、Ppp1r1b、Gprin3、Kcnab1、Pdyn、Gng7；PAG 腦 區 的 樞 紐 基 因 為Cd99l2、Smim17、Dzank1、Mthfd2l、Klhl1、Fuca2、Slc35g2、Pcyt1a、AW551984、Igsf1；mPFC 腦區的樞紐基因為Bmp3、Neurod6、Dlgap1、Adcy1、Olfm1、Pter、Tbr1、Id2、Sel1l3、Mical2。

討 論

慢性疼痛是臨床常見病，并常伴隨抑郁癥等情緒障礙。抑郁的病因和發病機制尚未完全清楚[11]。本研究用WGCNA 的方法分析Descalzi 等分享的RNA-seq 數據，對慢性疼痛導致抑郁動物的NAc、PAG 及mPFC 三個腦區的基因表達變化進行研究，共獲得9 個共表達模塊，將三個腦區組織作為研究性狀，計算與模塊的相關性，鑒定到三個腦區的特異性模塊，分析各自特異的基因表達情況，有助于探索三個腦區在該模型中的重要功能。

對特異性模塊內基因進行GO 功能注釋和KEGG 通路富集分析，發現這些基因主要富集在突觸結構、膜電位、離子通道、神經遞質等分子功能和生物學過程，反映了慢性疼痛引起了三個腦區的突觸結構和功能發生改變，是導致焦慮和抑郁的神經基礎。GO 分析發現在三個腦區都有小GTP 酶結合 (GO:0031267)、Ras GTP 酶結合 (GO:0017016) 兩個分子功能，其中Ras GTP 酶結合包含在小GTP酶結合功能中[12]。在突觸后致密區發現有多個小GTP 酶 (small GTPase) 亞家族蛋白的存在，參與調節樹突棘的形態，與突觸傳遞和突觸可塑性有關[13]。在抗抑郁治療的蛋白質組研究中發現，治療前有8 個Ras 相關蛋白表達有差異，治療6 周后有13 個Ras相關蛋白表達有差異，提示Ras 蛋白在抑郁中可能發揮重要作用[14]。同時，三個腦區之間存在神經纖維投射，共同進行疼痛和抑郁的調節。在急性和慢性疼痛期間mPFC 會發生神經遞質、基因表達、神經膠質細胞和神經炎癥的變化，從而導致其結構、活性和連接性發生改變，導致情緒決策受損[15]。向PAG 的皮質投射主要來源于mPFC，因而PAG 在疼痛中的功能受到mPFC 的調節。NAc 接受腹側被蓋區的多巴胺神經元投射和來自mPFC、海馬等腦區的谷氨酸能神經元投射，信息在NAc 進行整合。同時，NAc 經腹側蒼白球和黑質投射到丘腦，丘腦又投射到PFC、運動皮層及前扣帶回皮層，進行疼痛和情緒調節[16]。因此，慢性疼痛和抑郁癥在基因表達模式和神經可塑性機制上表現出共同的變化，并且彼此之間有很強的聯系。

不同腦區參與和執行不同的生理功能，因而基因表達具有組織特異性。本研究進一步通過Cytoscape 軟件分別對三個模塊進行分析，篩選出三個腦區各自的樞紐基因。NAc 的樞紐基因有Drd1a、Pdyn等。Drd1a編碼的Drd1受體，在獎賞通路中起重要作用，NAc 中D1 中型多棘神經元屬于直接通路，神經元活動性可被多巴胺增強，介導獎賞和正性情緒反應，而D2 中型多棘神經元屬于間接通路，神經元活動性可被多巴胺抑制，介導懲罰和負性情緒反應[17]。因而，D1 受體的激活可能通過增強獎賞和正性情緒來減弱疼痛，而D2 受體的激活可能會減少對疼痛的厭惡，二者的失衡會導致認知和情緒的異常[18]。有研究表明SNI 模型組小鼠經2 周抗抑郁藥地昔帕明治療后，NAc 中Drd1a表達增加，在抗抑郁中起重要作用[19]。Pdyn高表達于D1 中型多棘神經元，編碼的強啡肽是一種阿片肽，是κ 受體的特異性激動劑，可調節痛覺、抑郁、焦慮、成癮等多種生理功能[20]。PAG 的樞紐基因有Cd99l2、Pcyt1a等。Cd99l2是一種膜蛋白，被激活后介導白細胞通過血腦屏障進入神經元，引發神經系統免疫反應[21]，可能導致抑郁的發展。Pcyt1a編碼膽堿磷酸胞苷酰轉移酶A，在突觸蛋白磷脂合成方面起作用，而膜磷脂和突觸蛋白水平是突觸形成的間接指標[22]。mPFC 的樞紐基因有Neurod6、Dlgap1、Adcy1等。Neurod6調節神經分化與神經元存活，是認知相關的神經疾病的潛在治療靶點[23]。Dlgap1是一種突觸支架蛋白，位于谷氨酸能神經元的突觸后致密區(postsynaptic density, PSD) 處，是PSD95 組織的超分子蛋白復合物的組成部分，其突變會導致精神障礙和神經發育障礙[24]。Adcy1編碼的腺苷酸環化酶I 是一種神經特異性膜結合蛋白，催化關鍵的第二信使cAMP 的形成，Adcy1的表達與大腦中與可塑性相關的區域有關，包括海馬和大腦皮層[25]，其活性可能與抗抑郁藥的功效有關[26]。本研究中分析三個腦區得到了不同樞紐基因，提示在基因表達水平上，三個腦區在慢性疼痛導致的抑郁中有不同的作用。

目前對慢性疼痛的轉錄組研究主要集中在脊髓如背根神經節等部位的基因表達分析，近年來有研究開始關注慢性疼痛相關腦區的轉錄組研究。Mitsi 等[19]在Rgs9基因敲除小鼠構建SNI 模型，用RNA-seq 技術研究了NAc 內的RGS9-2 (regulator of G-protein signaling 9-2) 調節三環類藥物抗抑郁作用的機制，由于Rgs9基因敲除小鼠在SNI 建模后第5 周有顯著抗抑郁效果，進行轉錄組分析發現了Fezf2、Irs4、Myoc、Necab3等表達水平顯著變化的基因，但未發現與本研究相同的基因，這可能與取材時間有重要關系。Topham 等[27]在CD-1 小鼠構建SNI 模型，應用RNA-seq 技術在表觀遺傳水平研究了慢性疼痛發展過程中PFC 腦區的DNA 甲基化變化趨勢，GO 功能注釋發現在慢性疼痛的整個發展過程中很多基因涉及肌動蛋白和細胞骨架調節功能，表明外周損傷引起PFC 發生了持續的皮層可塑性和重塑，慢性疼痛導致PFC 功能逐漸異常。

抑郁有復雜的異質性，動物模型建模機制多樣以及各模型的抑郁癥狀表現不完全相同，使得不同抑郁模型的基因表達也不完全一致。SNI 模型與其他抑郁模型雖然存在大量不同的基因表達，但在文獻報道中也發現了部分表達一致的基因[28]，如SNI與慢性溫和應激模型 (chronic mild stress, CMS) 在mPFC 區有Hspb11、Cwc22、Nudt2基因表達上調，Ebp、Cbl基因表達下調；SNI 與慢性社會應激模型(chronic social stress, CSS) 在mPFC 區有Hspb11、Npas4、Xpc、Tipin基 因 表 達 上 調， 而Xbp1、Pold3、Sparc、Stx18、Mrpl38基因表達下調；SNI與CSS 模 型 在NAc 區 有Mpzl1、Tipin、Cep57、Mios、Ica1基因表達上調，Ythdf2、Cox17、RP23-86F8.3基因表達下調。

本研究將WGCNA 方法應用到慢性疼痛導致抑郁癥的相關腦區的基因表達研究，也存在不足之處。首先，分析得到的樞紐基因是否參與了慢性疼痛和抑郁的形成及其貢獻度，需在后續研究中進一步通過行為、分子生物學等實驗進行功能研究，如使用Real-time PCR 進行驗證后，使用基因修飾小鼠進行行為和功能驗證。其次，本研究采用的樣本量相對較小，這與目前數據庫中可用數據量有一定關系，需要在將來應用更多數據做進一步的深入分析。

綜上所述，本研究應用WGCNA 方法，分析了慢性疼痛導致抑郁的小鼠NAc、PAG 及mPFC三個腦區的基因共表達網絡，獲得了與腦區相關的基因模塊和樞紐基因，為此類抑郁癥的研究提供了潛在標志物，為該類疾病的發病機理研究提供了新思路。