夏博文，邢軍超，艾秋池，李瀚卿，徐美濤，侯天勇

陸軍軍醫(yī)大學(xué)第一附屬醫(yī)院骨科，重慶 400038

作為骨科常見病，下腰痛是威脅當(dāng)代人類身心健康的主要因素之一，椎間盤退行性變（IDD）則與很多下腰痛病例具有較高的相關(guān)性［1］。椎間盤是椎體間的連接結(jié)構(gòu)，其由位于中央的髓核（NP）、位于邊緣的纖維環(huán)（AF）以及與上下椎體毗鄰的軟骨終板（EP）構(gòu)成，它的主要功能是緩沖人自身體重及運動所產(chǎn)生的機械負荷［2］。然而，在退變狀態(tài)下，椎間盤將無法很好地滿足正常的生理需要［1］?，F(xiàn)階段，保守和手術(shù)方法均可用于IDD的臨床治療，前者包括對癥的非甾體抗炎藥（NSAIDs）、物理療法和行為調(diào)整等，后者則主要是椎間盤摘除及椎體融合術(shù)［3-5］，但手術(shù)病人將承受數(shù)月的受限腰部活動和可能發(fā)生的鄰椎病等并發(fā)癥［6］。因此，為了改進臨床治療方法和改善病人預(yù)后，尋找IDD的關(guān)鍵影響因素并發(fā)展對因治療方法便成了此領(lǐng)域的主要研究方向。

轉(zhuǎn)錄組測序（RNA-seq）技術(shù)可用來獲取某物種細胞或組織的完整轉(zhuǎn)錄組數(shù)據(jù)，從而揭示該物種在不同生長環(huán)境下的基因表達變化［7］。既往已有學(xué)者使用此方法分析IDD進展過程中的差異基因和通路，但由于臨床正常椎間盤樣本難以獲取，他們將椎間盤突出患者樣本和椎體滑脫（患者樣本進行對比，具有一定誤差［8］。同時，動物模型也被應(yīng)用于相似的研究中，然而由于人與動物間不同的分子表型和椎間盤結(jié)構(gòu)，其結(jié)論仍有爭議［9］。

本課題組通過RNA-seq技術(shù)，首先對臨床手術(shù)采集到的椎間盤樣本進行測序，其次根據(jù)其有無IDD征象分為2組，對比所得的測序數(shù)據(jù)，得到差異基因，并根據(jù)它們在通常生理過程中的不同作用，分別闡述它們在實驗結(jié)果中的表達特征和可能意義。我們隨后利用GO和KEGG數(shù)據(jù)庫將上述差異基因富集到相應(yīng)通路中，結(jié)合既往相關(guān)研究結(jié)果，分析這些通路在臨床工作或后續(xù)研究中可能起的作用。

1 材料和方法

1.1 樣本收集

椎間盤髓核樣本收集自臨床進行手術(shù)的椎間盤突出/脫出病人（n=4，年齡53～68歲，2男1女，其中兩份樣本來自于同一病人的不同節(jié)段）和臨床因為治療需要而切除正常間盤的患者，如需要行椎間融合手術(shù)的胸椎管狹窄癥、脊柱側(cè)彎的患者（n=3，年齡20～48 歲，2 男1女）。上述受試者的IDD進展程度均經(jīng)過MRI檢查，并由改良Pfirrmann退變標準評估分級［10］，其中IDD組患者的病變節(jié)段經(jīng)評估分別為5、6、7、8級，非IDD組則均不超過2級（圖1）。

圖1 部分入組患者的MRI檢查影像Fig.1 MRI images in selected cases.A,B:T2-weighted MR images confirming non-degenerated L5/S1 IVD in a case of spondylolisthesis case and T10/T11 IVD in a case of thoracic spinal stenosis case.C:Obvious degeneration of L4/L5 IVD in a case of intervertebral disc herniation.Arrows indicate the segments for surgery.

IDD組的排除標準包括（1）年齡<40歲或>70歲；（2）合并峽部骨折的腰椎滑脫癥；（3）合并其它系統(tǒng)疾病，如嚴重呼吸系統(tǒng)疾病，心腦血管疾病，消化系統(tǒng)疾病，內(nèi)分泌系統(tǒng)疾病，泌尿系統(tǒng)疾病，凝血功能障礙等；（4）精神病患者，有吸毒、酗酒不良嗜好者，懷孕婦女；（5）同時參加其他研究者或剛結(jié)束其他臨床試驗研究者。非IDD組的排除標準則為前者的（3）、（4）、（5）兩點。本課題研究通過中國重慶市陸軍軍醫(yī)大學(xué)西南醫(yī)院倫理委員會審批，受試者均經(jīng)過知情同意（批準編碼：KY2019141）。

1.2 RNA提取及評估

所有樣本的RNA均按照制造商指南使用Trizol試劑（Invitrogen）提取，并立即隔離儲存于-80 ℃環(huán)境。隨后RNA樣本按說明書指導(dǎo)經(jīng)Direct-zol RNA miniPrep Kit（ZYMO research）試劑盒分離，RNA 濃度及260/280、260/230吸光度比值則使用NanodropOne分光光度計（Thermo）進行測量［11］。

1.3 高通量轉(zhuǎn)錄組測序

本檢測項目采用EASY RNA-seq方法完成7個樣本的建庫，經(jīng)質(zhì)控檢測后，采用Illumina平臺測序［11］。樣本測序數(shù)據(jù)質(zhì)量的評估標準包括：樣本原始測序數(shù)據(jù)的Q20百分比需超過90%，樣本原始測序數(shù)據(jù)過濾比例大于80%。接下來使用HISAT2比對樣本過濾后的測序讀長和參考基因組［12］，再使用Bowtie2比對樣本過濾后的測序讀長和參考基因集［13］，人類基因的參考序列文件和基因組注釋文件均來源GENCODE數(shù)據(jù)庫［14］。通過上述質(zhì)量檢測的樣本測序數(shù)據(jù)則進行接下來的生物信息學(xué)分析流程。

1.4 生物信息學(xué)分析

首先使用FeatureCounts進行基因原始表達量的統(tǒng)計［15］，為了消除樣本測序數(shù)據(jù)量、基因序列長度等差異引起的基因表達量統(tǒng)計偏差，我們采用了FPKM（Fragments Per Kilobase of exon model per Million mapped fragments）方法對樣本基因表達量進行標準化［16］，并采用小提琴圖（Violin Plot）繪制樣本基因表達從而標準化數(shù)據(jù)的分布情況。接著使用R 工具包DESeq2分析組間差異表達基因，在本研究中，顯著表達差異篩選標準設(shè)定為：P值經(jīng)過多重校驗校正后的值（padj）≤0.05［17］。

我們隨后繪制熱圖對上述實驗流程所獲得的差異基因進行可視化并利用K-MEANS聚類方法進行分類，從而更好地研究差異基因的不同表達模式［18］。最后我們分別利用GO和KEGG數(shù)據(jù)庫對上述結(jié)果進行富集注釋，即通過與整個基因組背景相比，分析差異表達基因列表顯著富集的相應(yīng)數(shù)據(jù)庫條目。

1.5 qRT-PCR驗證

我們使用GoTaq?qPCR Master Mix 和ABI ViiATM7 Real-Time PCR System驗證了組件差異基因的表達。qRT-PCR引物由Primer 5.0軟件設(shè)計，引物詳情見表1。管家基因3-磷酸甘油醛脫氫酶(GAPDH)用作內(nèi)參對照。完整的熱循環(huán)體系包括95 ℃3 min、95 ℃15 s以及60 ℃和72 ℃各30 s，共計40循環(huán)，熔解階段由95 ℃30 s和60 ℃1 min組成，運用2-ΔΔCt法計算不同組別目標mRNA間的倍數(shù)變化［19］，數(shù)據(jù)以差異倍數(shù)的log2值（log2FoldChange）表示，使用R-4.0.4進行分析，采用Wilcoxon符號秩檢驗進行統(tǒng)計學(xué)處理。設(shè)定P<0.05為差異有計學(xué)意義。

表1 引物信息Tab.1 The primer information

2 結(jié)果

2.1 轉(zhuǎn)錄組測序及建庫

本次研究分別從IDD及非IDD組的RNA文庫測得總計84456816和44816454條原始讀長，經(jīng)過篩選，分別剩下78333214和42085446條可用讀長，這些可用讀長則進一步與人類基因組數(shù)據(jù)進行對比從而進行下一步的數(shù)據(jù)篩選（表2）。

表2 測序數(shù)據(jù)總結(jié)Tab.2 Summary of sequencing data

2.2 組間相關(guān)性研究及差異基因的獲取

差異基因的顯著性指標設(shè)置為padj≤0.05，將符合此標準的差異基因標注在組間散點圖中從而更好地觀察它們的散布及表達特征。MA圖則用來展示顯著差異基因表達水平和倍數(shù)變化的散布模式?；鹕綀D展示了兩組間最具顯著性的差異基因（|log2FoldChange|≥1且padj≤0.05）。按上述標準，相對于非IDD樣本，共得到了210個上調(diào)表達、302個下調(diào)表達的顯著差異基因（圖2）。

圖2 組間差異基因的篩選Fig.2 Identification of intergroup DEGs.A:Every point in the figure stands for a gene,X and Y axes separately represents the gene expression level of different groups.Gray dots represent genes with insignificant expression difference,orange dots represent genes which are significantly up-regulated in IDD group compared with non-IDD group,and green dots represent genes that are significantly down-regulated.B:MAdiagram.Each point in the figure represents a single gene,and the X axis represents the average expression level.The Y-axis represents log2FoldChange (IDD/non-IDD expression).C:Volcano plots.Each point in the figure represents a gene,X axis:log2FoldChange,Y axis:-log10(padj).The red dots:DEGs with whose|log2FoldChange|≥1 and padj≤0.05;the green dots:DEGs with whose padj>0.05 and|log2FoldChange|≥1;the grey dots are the non-DEGs.The red-dot-genes were taken into the further analysis.

2.3 差異基因的富集分析

根據(jù)分組信息計算不同基因的組內(nèi)平均表達量，然后用Z評分對其進行標準化。隨后對上述標準化平均表達量進行層次聚類，結(jié)果以熱圖的形式展示并基于上述結(jié)果進行K-MEANS聚類分析。正常組樣本和異常組樣本的基因表達趨勢不同，而組內(nèi)的基因表達趨勢較一致。另外，表達上調(diào)的差異基因雖然數(shù)量少于下調(diào)基因，但其各基因間的歐氏距離相對較大，上調(diào)基因具有較高的離散程度，相對難以富集到一致的通路上（圖3）。

圖3 差異基因的富集分析Fig.3 Enrichment analysis of the DEGs in IDD.A:Differential gene cluster analysis diagram.The column represents the expression profile of DEGs,and the row represents the genetic expression quantity among different samples,the left 3 are the samples from non-IDD group,and the right 4 are IDD samples.B:K-means cluster diagram of intergroup DEG set,which also shows the Euclidean distance between genes.

2.4 差異基因的功能富集

分別展示了GO數(shù)據(jù)庫生物學(xué)過程、細胞組分、分子功能3個模塊以及KEGG數(shù)據(jù)庫的差異基因功能富集結(jié)果（圖4）。

圖4 差異基因的GO和KEGG聚類分析Fig.4 GO and KEGG enrichment analysis of the DEGs.A-C:The top 20 terms summarized from the GO biological process,cellular components and molecular function databases.D:Result of KEGG enrichment.

2.5 使用qRT-PCR驗證轉(zhuǎn)錄組測序結(jié)果

根據(jù)log2FoldChange絕對值的大小、預(yù)期研究價值和既往文獻支持選取了13個差異基因作為qRT-PCR的驗證目標，分別是FN1、CXCL8、CXCL1、IGFBP6、TNFAIP6、TIMP1、TIMP2、PF4、COL1A1、COL3A1、POSTN、CHI3L2和MMP3，它們在兩組間表達具有差異顯著性（P<0.001），其中TIMP1、TIMP2和CXCL8的qRT-PCR驗證結(jié)果經(jīng)2-ΔΔCt法計算樣本組間差異無顯著性（Log2FoldChange≤1），其余10個基因具有差異顯著性（Log2FoldChange>1）且和RNA-seq相似的表達趨勢（圖5）。

圖5 轉(zhuǎn)錄組測序結(jié)果的qRT-PCR驗證Fig.5 Verification of RNA-Seq data by qRT-PCR.Ten of the 13 DEGs show significant intergroup fold change(Log2FoldChange>1),whose expression trends are consistent with the results of RNA-seq.

3 討論

炎癥反應(yīng)一直是影響包括IDD在內(nèi)的多種退行性疾病發(fā)生發(fā)展的重要因素，在本次研究中多種炎癥因子也在組間具有顯著性差異表達。白介素在機體炎癥介導(dǎo)和調(diào)控方面具有重要地位，在本次的轉(zhuǎn)錄組差異基因中，注意到包括IL32、IL1B、IL2RA、IL7R、IL1RN、IL36RN、IL22RA2、IL1A、IL36G在內(nèi)的多個白介素相關(guān)編碼基因?qū)Ρ日Ｗ甸g盤組織表達均呈下調(diào)趨勢，其中既包括具有促進炎癥反應(yīng)的因子，也有IL36RN等炎癥抑制劑。這樣的結(jié)果似乎與既往類似的研究結(jié)果不符，一般認為IDD 會伴有大部分炎癥介質(zhì)的表達上調(diào)［20］。然而IDD是一個長期的過程，在病程較長的椎間盤突出病例中（距首次癥狀出現(xiàn)超過半年），炎癥反應(yīng)反而會因退變組織的部分萎縮和降解逐漸減弱。與之類似，C-X-C模體趨化因子配體（CXCL）與C-C模體趨化因子配體（CCL）家族基因，包括CXCL10/9/2/3/1/8 和CCL5/19/22作為人體重要的趨化因子成分之一在IDD組椎間盤樣本內(nèi)表達相對正常組也均呈下調(diào)趨勢。椎間盤是人體最大的無血管組織，其與血源性的免疫因子接觸十分有限，因此，在生理狀態(tài)下正常的免疫炎癥反應(yīng)很難作用于椎間盤組織。然而在椎間盤突出發(fā)生時，髓核通過破裂的纖維環(huán)與椎管接觸，并且局部血運在生血管相關(guān)基因的作用下增加，增大了炎癥趨化因子與退變的椎間盤組織作用的可能性。通過炎癥因子的作用，突出的椎間盤組織會導(dǎo)致疼痛反應(yīng)，在少數(shù)情況下，其也可能被機體吸收［21］。在本次建立的RNA-seq文庫中，多種炎癥因子的下調(diào)表達可能反映了退變晚期，突出椎間盤無法被機體自我吸收時的基因表達情況。

細胞外基質(zhì)組分是維持椎間盤髓核細胞正常生態(tài)的重要成分，其對維持微環(huán)境的含水量及正常代謝環(huán)境具有重要作用。膠原蛋白是椎間盤細胞外基質(zhì)的主要成分之一，已有多個研究提示其不同類型膠原蛋白的含量變化與IDD進展情況具有明顯相關(guān)性，尤其是Ⅱ型膠原蛋白含量的減少和Ⅰ型膠原蛋白含量的升高被認為是髓核細胞退變的重要標志之一［22］。在本次建庫數(shù)據(jù)中，膠原蛋白編碼基因COL1A1、COL4A1、COL3A1、COL4A2、COL18A1和COL12A1均是顯著差異基因且對比正常組呈表達升高趨勢，而COL2A1雖然未達顯著性差異標準，但對比表達下調(diào)，與預(yù)期結(jié)果基本一致，其中部分膠原蛋白成分的可能角色在過去并未被報導(dǎo)。其他細胞外基質(zhì)組分，包括細胞角蛋白（KRT）等成分，作為正常髓核細胞的表型標志，也呈下調(diào)表達［23］?；|(zhì)金屬蛋白酶（MMP）具有強烈的基質(zhì)降解作用，一直被認為是IDD的重要引導(dǎo)因素，其存在已被證實可能加重包括骨關(guān)節(jié)炎在內(nèi)的多種退行性疾?。?4］。本次測序證實了MMP3、MMP2在退變椎間盤組織中的上調(diào)表達，以及其抑制因子TIMP1/2/3則反射性地表達增強。除上述外，還在差異基因文庫中找到了其他可能重要的細胞外基質(zhì)成分編碼基因。幾丁質(zhì)酶3樣蛋白2（CHI3L2）缺乏幾丁質(zhì)酶活性，并已被證實在機體對組織重組、損傷和炎癥反應(yīng)方面具有關(guān)鍵地位，其在骨關(guān)節(jié)炎發(fā)生進展中調(diào)控軟骨再生和Ⅱ型膠原纖維蛋白穩(wěn)態(tài)的作用尤其值得重視［25］。近期有學(xué)者指出，它的同源基因CHI3L1可能通過AKT3信號通路在IDD進展過程中起保護性作用，但其中的具體調(diào)控機制尚不明確［26］。骨膜蛋白（POSTN）是一類廣泛存在于骨骼、皮膚、韌帶等人類結(jié)締組織中的細胞外基質(zhì)分泌蛋白［27］，體外實驗已證實其可與MMP13基質(zhì)金屬蛋白酶協(xié)同作用從而參與骨關(guān)節(jié)炎進程中的軟骨基質(zhì)降解過程［28］。另有研究者通過免疫學(xué)和蛋白組學(xué)實驗發(fā)現(xiàn)POSTN在退變椎間盤細胞中的異常激活，但未進行進一步的探究［29］。以上2種細胞外基質(zhì)構(gòu)成組分基因在本次差異對比中上調(diào)表達且展現(xiàn)了很高的差異度--其log2FoldChange分別為5.566691773和4.28670097，位于所有差異基因的前10。

眾多的生長因子廣泛參與人機體的各項生物過程中，椎間盤組織細胞的代謝調(diào)控也不例外。既往已有文獻指出，包括胰島素樣生長因子（IGF）和轉(zhuǎn)化生長因子（TGF）在內(nèi)的多個生長因子家族基因變異可能與IDD有關(guān)［30］。本次測序發(fā)現(xiàn)，退變組織中的胰島素樣生長因子結(jié)合蛋白（IGFBP)4/5/6/7對比未退變組織均有不同程度地表達上調(diào)，然而與此有所矛盾的是既往研究結(jié)果中IGFBP5被認為主要在椎間盤組織中發(fā)揮促細胞增殖、抑制細胞凋亡的作用，且往往在IDD進展過程中表達逐漸下調(diào)，而全外顯子測序結(jié)果則猜測IGFBP6的變異與IDD有較高的相關(guān)性［31］，或許IGF家族及其受體間存在更復(fù)雜的信號通路促進或抑制IDD的進展。TGFβ是人體生長發(fā)育過程中最活躍、作用最為廣泛的信號通路中，各項體內(nèi)體外實驗已證實相當(dāng)數(shù)量的編碼基因通過與該通路復(fù)雜的交互作用從而影響椎間盤的微環(huán)境變化［32］。

隨后對上述富集到的所有差異基因進行了GO和KEGG數(shù)據(jù)庫功能富集分析。GO數(shù)據(jù)庫可以從3個方面--生物學(xué)過程，細胞組分和分子功能--分析兩樣本組間差異基因，從而幫助了解IDD在不同生物學(xué)層面上的影響因素，從而更加具體地制定下一步研究計劃或提出針對性臨床診療方案。本次主要分析了在各個模塊富集基因占比位于前20的通路。在GO數(shù)據(jù)庫生物過程模塊的富集結(jié)果中，約一半的通路與上皮細胞發(fā)育及角化相關(guān)，余下的通路則主要包括炎癥趨化和細胞外基質(zhì)合成，綜合各通路中富集到的具體差異基因，前者成為差異通路可能與椎間盤中豐富的纖維成分和NP細胞的類上皮細胞表型在IDD中的異常改變有關(guān)［33］。細胞組分方面，細胞外基質(zhì)相關(guān)組分無論是在富集度還是基因差異顯著度方面都位于前列，而其余通路則包含了各類常見細胞器。綜上，這兩個模塊的富集結(jié)果與既往研究結(jié)果基本一致，并再次強調(diào)了細胞外基質(zhì)、胞外膠原蛋白成分和炎癥因子在IDD發(fā)生進展過程中所起的重要作用，后續(xù)研究仍可圍繞這些成分展開。近來，分子功能越來越成為多種腫瘤、遺傳病的研究重點，因此，本研究期望GO數(shù)據(jù)庫的分子功能富集分析模塊能帶來更細致的答案。在所有富集到的通路中，生長因子結(jié)合（GO:0019838）和細胞外基質(zhì)結(jié)構(gòu)構(gòu)成（GO:0005201）的顯著度遠高于其他通路。此前我們已描述了生長因子家族在IDD發(fā)生、發(fā)展甚至在抑制其發(fā)展方面的復(fù)雜作用。另外，基于椎間盤細胞質(zhì)量較少，富含水分和膠原蛋白的環(huán)境，細胞外基質(zhì)對維持NP細胞穩(wěn)態(tài)、負荷外界機械力不可或缺［33］。其余富集到的差異基因通路還包含了多種趨化因子、整合素、肝素與相應(yīng)受體的結(jié)合過程等。KEGG數(shù)據(jù)庫更注重分析基因產(chǎn)物在細胞中的代謝途徑以及這些基因產(chǎn)物功能，在本次的分析結(jié)果中，不僅IL-17、TNF、PI3K-Akt等對IDD具有調(diào)控作用的經(jīng)典信號通路顯示顯著性富集作用，而AGE-RAGE和部分糖尿病相關(guān)通路也有一定差異性。令人意外的是KEGG富集通路中差異顯著性最高的2條分別為阿米巴病（hsa05146，Amoebiasis）和病毒蛋白與細胞因子及其受體相互作用（hsa04061），此外，軍團菌（hsa05134，Legionellosis）和百日咳（hsa05133，Pertussis）通路也被差異基因顯著性富集，暗示IDD與外界病原體感染有關(guān)。然而，IDD作為一種退行性疾病，學(xué)界公認其發(fā)生發(fā)展與年齡增長、機體功能衰退及外界負荷增大有關(guān)，一般不會認為病原體對IDD的發(fā)生進展會起到首要作用。既往有少數(shù)研究指出病原體隱形感染可能會推動IDD進展，但缺乏更多體內(nèi)實驗證實，本次的發(fā)現(xiàn)可能為類似觀點提供了新論據(jù)［34］。

本研究有如下的缺陷：（1）樣本量較小，無法很好地更好地設(shè)置生物學(xué)重復(fù)并排除個體間差異。另外，入組患者間，包括年齡、性別、不同的椎間盤節(jié)段等也沒有很好地統(tǒng)一。造成樣本量不足的主要原因是臨床上正常椎間盤樣本的獲取困難--隨著診療方案和手術(shù)技術(shù)的進步，非IDD手術(shù)患者都越來越少地選擇椎間盤摘除術(shù)。下一步本課題組可能會使用動物模型來部分彌補該缺陷；（2）RNA-seq技術(shù)盡管可以精確快速地獲得某一物種特定細胞或組織在某一狀態(tài)下的所有轉(zhuǎn)錄本的集合，從而深入研究mRNA和非編碼RNA在IDD進展過程中的變化情況，但其分析僅局限于某一特點時間點，再加上很難確保入組病人處于完全相同的IDD進展階段，因此可能會帶來額外的誤差；（3）手術(shù)及手術(shù)取樣時往往是切除整塊髓核組織進行RNA提取及后續(xù)分析，無法細致地分析髓核內(nèi)不同細胞群的IDD發(fā)生進展因素差異，下一步考慮采用單細胞測序技術(shù)，具體區(qū)分椎間盤內(nèi)的各細胞群組，以研究各細胞系對IDD的獨特影響；（4）本次研究沒有對退變程度不同的樣本進行分類差異分析，無法得出IDD進展不同階段內(nèi)各自最關(guān)鍵的影響因素，繼而尋找IDD進展到不同程度時最佳的治療方案，這與有限的樣本含量和RNA-seq的時限性特點有關(guān)。下一步研究本課題組擬加大樣本量，并根據(jù)改良Pfirrmann分級方法和術(shù)中所見病人的突出/脫出程度對入組病人進行分類，以求減少誤差，獲取更加精確的結(jié)論。