軟骨細胞單細胞轉(zhuǎn)錄組數(shù)據(jù)的再挖掘揭示人骨關(guān)節(jié)炎病程中轉(zhuǎn)錄調(diào)控因子和代謝的時序變化

施劍明,杜桂花

(1.景德鎮(zhèn)市第一人民醫(yī)院創(chuàng)傷骨科,江西 景德鎮(zhèn) 333000;2.南昌大學(xué)公共衛(wèi)生學(xué)院勞動衛(wèi)生與毒理學(xué)教研室,江西省預(yù)防醫(yī)學(xué)重點實驗室,南昌 330006)

人骨關(guān)節(jié)炎(human osteoarthritis,hOA)是導(dǎo)致老年人運動能力喪失的主要原因之一,給社會和家庭造成巨大的經(jīng)濟損失[1-2]。誘導(dǎo)hOA發(fā)生的因素眾多且關(guān)系復(fù)雜,有研究[3]表明,代謝、遺傳、衰老和損傷等都與hOA的發(fā)生發(fā)展相關(guān)。目前,臨床對hOA的治療主要以非類固醇類抗炎治療和物理治療為主,但都無法逆轉(zhuǎn)疾病進程,人工關(guān)節(jié)置換通常是患者的最終選擇[4]。探究hOA的發(fā)病及進展機制,可為逆轉(zhuǎn)病情進展治療方案的開發(fā)提供方向和數(shù)據(jù)支持,是目前的研究熱點。

臨床上,根據(jù)國際軟骨修復(fù)協(xié)會(International Cartilage Repair Society,ICRS)軟骨損傷分級系統(tǒng)和國際骨關(guān)節(jié)炎研究協(xié)會(Osteoarthritis Research Society International,OARSI)評分標(biāo)準(zhǔn),骨關(guān)節(jié)炎處軟骨細胞受損程度可分為由輕微到嚴(yán)重的4個等級:正常無損傷或軟骨表面完整處為0級(Stage 0,S0);表淺的鈍性損傷或軟骨表面不平整處為I級(Stage I,S1);軟骨損傷小于軟骨厚度一半或軟骨表面不連續(xù)處為Ⅱ級(Stage Ⅱ,S2);軟骨損傷大于軟骨厚度一半但未達軟骨下骨或軟骨出現(xiàn)垂直裂變處為Ⅲ級(Stage Ⅲ,S3);軟骨損傷可見軟骨下骨暴露或出現(xiàn)軟骨侵蝕處為Ⅳ級(Stage Ⅳ,S4)。目前,盡管大量研究[5-7]報道了hOA的病理表現(xiàn)和潛在分子機制,但并未對其病情進程中的分子動態(tài)變化進行探索。因此,對hOA進展過程中的調(diào)控網(wǎng)絡(luò)進行研究可為臨床阻止或逆轉(zhuǎn)hOA進程提供理論基礎(chǔ)。

單細胞轉(zhuǎn)錄組測序(single-cell RNA sequencing,scRNA-seq)利用有限數(shù)量的樣本在單細胞水平上獲得轉(zhuǎn)錄組信息,可為發(fā)育生物學(xué)的研究和生理病理機制的探索提供重要信息[8-9]。JI等[10]通過scRNA-seq分析了hOA患者軟骨細胞的內(nèi)部異質(zhì)性及分子特征,但并未從hOA病程分級角度上分析疾病進展過程中基因表達及其生物學(xué)過程改變。而本研究通過對人膝關(guān)節(jié)軟骨細胞的單細胞測序數(shù)據(jù)進行再挖掘,從單細胞水平上對hOA各發(fā)展階段中基因及生物學(xué)過程變化進行連續(xù)性的研究,以期為臨床開展hOA的階段針對性治療和新療法的開發(fā)建立基礎(chǔ)。

1 數(shù)據(jù)和方法

1.1 原始數(shù)據(jù)來源

人軟骨細胞單細胞測序數(shù)據(jù)集(序列號:GSE104782)[10]下載于GEO數(shù)據(jù)庫。該數(shù)據(jù)集中軟骨細胞樣本來源于10例膝關(guān)節(jié)置換手術(shù)患者的膝關(guān)節(jié),依據(jù)每位患者關(guān)節(jié)損傷情況分別采集S0—S4分級的軟骨細胞進行單細胞轉(zhuǎn)錄組建庫及測序。

10例患者中,男3例、女7例,平均年齡(66.50±1.86)歲,體重(66.90±2.44)kg,美國特種外科醫(yī)院膝關(guān)節(jié)評分(hospital for special surgery knee score,HSS)(60.20±2.56)分,具體見表1。

1.2 數(shù)據(jù)挖掘方法

基于R軟件對單細測序數(shù)據(jù)集進行再挖掘。

1.2.1 質(zhì)量控制

進行數(shù)據(jù)分析之前,利用單細胞檢測mRNA分子總數(shù)(nCount_RNA)和基因數(shù)量(nFeature_RNA)剔除測序深度不足或因細胞粘連導(dǎo)致測序質(zhì)量不佳的單細胞測序數(shù)據(jù)。據(jù)圖用于后續(xù)分析的準(zhǔn)入標(biāo)準(zhǔn)為:“2000

1.2.2 數(shù)據(jù)分析和可視化

基于R 3.6軟件,通過軟件包Seurat 3.0和Monocle 2進行數(shù)據(jù)分析和可視化的實現(xiàn),其中在Seurat 3.0中進行測序數(shù)據(jù)的初始標(biāo)準(zhǔn)化和歸一化,在Monocle 2中進行細胞擬時軌跡分析。通過軟件包CellCycleSoring進行細胞周期檢測,通過DAVID Bioinformatics Resources 6.8(https://david.ncifcrf.gov/)進行差異基因功能(gene ontology,GO)分析,通過PANTHER分類系統(tǒng)(http://pantherdb.org/)尋找差異基因中編碼轉(zhuǎn)錄調(diào)控因子的基因。

使用GraphPad Prism 8.0、R 3.6和Adobe Illustrator CC 2018等軟件進行圖表制作。

1.3 統(tǒng)計學(xué)方法

2 結(jié)果

2.1 單細胞測序數(shù)據(jù)概況

對原始測序數(shù)據(jù)進行質(zhì)控分析后,將數(shù)據(jù)過濾指標(biāo)設(shè)定為“2000

對篩選后的細胞進行標(biāo)準(zhǔn)化和歸一化分析后,進行均勻流形近似和投影(uniform manifold approximation projection,UMAP)降維分析(圖1A)。在此基礎(chǔ)上進一步驗證軟骨細胞已知標(biāo)志基因(如ACAN、COL2A1、ITM2A和VEGFA)的表達(圖1B),結(jié)果明確了篩選后的細胞均為軟骨細胞。雖然來自于10例患者的軟骨細胞并未各自成簇,但樣本間仍有異質(zhì)性(圖1C)。為檢測單細胞處理的批次差異效應(yīng),通過修正Pearson相關(guān)系數(shù)評定各樣本間相似程度[11],結(jié)果顯示Pearson相關(guān)系數(shù)范圍在0.899～0.983,樣本間雖存在異質(zhì)性,但一致性程度較高,可合并用于后續(xù)分析(圖1D)。

A:軟骨細胞scRNA-seq文庫UMAP;B:軟骨細胞標(biāo)志蛋白UMAP;C:10例hOA患者軟骨細胞scRNA-seq文庫UMAP;D:10例患者文庫間Pearson相關(guān)系數(shù)熱圖。

2.2 hOA發(fā)展進程中軟骨細胞的擬時軌跡分析

利用Monocle 2從單細胞水平上擬合hOA進展軌跡,并對隨軌跡變化的相關(guān)基因進行分析,以從單細胞水平上分析hOA發(fā)展的調(diào)控網(wǎng)絡(luò)。

各分級軟骨細胞UMAP結(jié)果顯示,雖然同一分級內(nèi)的細胞聚集成簇且S0至S4排列具有一定方向性,但并非所有細胞都按該規(guī)律排布(圖2A)。單細胞測序數(shù)據(jù)來自10例患者,雖然樣本間相關(guān)性分析提示樣本間的一致性較好,但不同分級軟骨細胞的采集是依據(jù)關(guān)節(jié)軟骨病損程度進行的。為明確各樣本分級,有學(xué)者[12-13]檢測了與骨關(guān)節(jié)炎進展呈正相關(guān)的細胞外基質(zhì)降解相關(guān)基因,結(jié)果顯示,ADAMTS12、ADAMTS4和ADAMTS16在hOA患者軟骨細胞中的表達升高,且S4軟骨細胞中以上基因表達水平均最高(圖2B),這表明從hOA進展的角度看各樣本一致性仍較好,可基于此進行擬時軌跡分析(pseudo-time trajectory analysis)。

A:不同分級(S0-S4)軟骨細胞UMAP;B:細胞外基質(zhì)降解相關(guān)標(biāo)志基因在不同分級軟骨細胞中的表達水平,其中點的大小表示基因表達的細胞百分比(Percent Expressed),灰度梯度表示基因均化表達量(Average Scaled Expression);C:各分級軟骨細胞沿著擬時軌跡軸排列;D:沿擬時軌跡軸變化的差異基因及其表達熱圖;E:各組差異基因GO分析部分條目;F:各組差異基因KEGG通路分析;G:根據(jù)各組差異基因及其表達峰值總結(jié)hOA各發(fā)展階段生物學(xué)過程。

利用各分級細胞前1500個差異基因(基于P值)將細胞進行排序,除個別S4的細胞排列在軌跡一端起點外,其余細胞根據(jù)其分級很好地依擬時軌跡軸排列,表明擬時軌跡分析模型構(gòu)建較成功(圖2C)。在此基礎(chǔ)上,進一步將沿著擬時軌跡軸變化的差異基因(P<0.01)及其表達量擬合成熱圖(圖2D),結(jié)果顯示差異基因可分為4個組,組1的基因表達在擬時軌跡軸最開始處表達最高,組2的基因在擬時軌跡軸前中段表達最高,組3基因在擬時軌跡軸后中段表達最高,而組4基因在擬時軌跡軸后段表達最高。

對各組差異基因進行GO分析和KEGG分析,細胞對TNF的應(yīng)答/TNF信號通路、對胰島素的應(yīng)答/MAPK信號通路、MAPK級聯(lián)反應(yīng)/PI3K-Akt信號通路分別出現(xiàn)在組1,組2和組3,而在組4中則出現(xiàn)細胞黏附和蛋白的消化、吸收過程(圖2E—F)。差異基因的表達峰值呈現(xiàn)從組1至組4的前后排列,因此基于差異基因可歸納總結(jié)出hOA發(fā)展4個階段的相關(guān)生物學(xué)過程(圖2G):在hOA發(fā)展的早期階段,其生物學(xué)過程主要涉及細胞因子的產(chǎn)生和一氧化氮的生物合成等;而hOA發(fā)展的晚期階段,關(guān)節(jié)軟骨組織中出現(xiàn)血管生成、膠原纖維形成和基質(zhì)破壞等;炎癥反應(yīng)則貫穿hOA進展的整個過程。

2.3 調(diào)控hOA發(fā)展的轉(zhuǎn)錄調(diào)控因子

包括RUNX2、MAFF和ATF3等在內(nèi)轉(zhuǎn)錄調(diào)控因子的在hOA的發(fā)生和發(fā)展過程中起重要作用[14-16]。為進一步從轉(zhuǎn)錄調(diào)控水平解析hOA進展中的關(guān)鍵因子,本研究利用PANTHER分類系統(tǒng)對隨著hOA軟骨細胞擬時軌跡變化的各組基因進行了分析(圖3)。已知的對hOA存在關(guān)鍵調(diào)控作用的轉(zhuǎn)錄因子ATF3、MAFF和RUNX2分別出現(xiàn)在組1,組2和組4。而FOXO3、KLF10、KMT2C和ZMYM2在hOA的發(fā)展初期表達最高(圖3A),轉(zhuǎn)錄調(diào)控因子HES1、ID1、ID3、JUN、KLF2和MEF2C表達則隨hOA擬時軌跡先下調(diào),隨后又在擬時軌跡后期表達上調(diào)(圖3B);轉(zhuǎn)錄調(diào)控因子ETS2、FOXA3、JUND、NFATC1和SOX8在hOA擬時軌跡的前中段表達最高(圖3C),而CITED4、EGR2、ETV5、FHL2和SOX11則在擬時軌跡的中后段表達最高(圖3D);轉(zhuǎn)錄調(diào)控因子AEBP1、CKLF和STAT1在hOA擬時軌跡后期表達上調(diào),這可能與hOA晚期軟骨組織自我修復(fù)可能性的徹底喪失有關(guān)(圖3E)。

A—B:組1;C:組2;D:組3;E:組4。圖3 hOA軟骨細胞中各組關(guān)鍵轉(zhuǎn)錄調(diào)控因子沿擬時軌跡軸的表達模式

2.4 hOA發(fā)展過程中軟骨細胞的能量代謝改變

對沿hOA進展擬時軌跡變化的軟骨細胞差異基因進行KEGG分析,結(jié)果顯示氧化磷酸化相關(guān)基因出現(xiàn)在組4中,且相關(guān)基因的表達隨hOA進展升高(圖4A)。關(guān)節(jié)軟骨是一種無血管組織,這意味著軟骨細胞對營養(yǎng)物質(zhì)的攝取及其代謝產(chǎn)物的排出主要依賴于軟骨表面的滲透,而軟骨細胞主要通過糖酵解過程獲得能量[17]。細胞能量代謝與細胞周期相關(guān),增殖活躍的細胞(如腫瘤細胞)更傾向于利用糖酵解獲取能量[18]。為排除hOA軟骨細胞中能量代謝改變是由于細胞增殖變化導(dǎo)致的,本研究利用CellCycleSoring對軟骨細胞的細胞周期進行了分析。結(jié)果表明,隨著hOA進展,處于增殖狀態(tài)(S/G2/M期)的軟骨細胞百分比明顯減少(圖4B),排除了hOA進展后期軟骨細胞氧化磷酸化水平增加是由細胞增殖活力改變而引起的可能性?？紤]到個別S4細胞同大部分S0細胞一樣排列于擬時軌跡軸一端起點,將擬時軌跡分析差異基因中與氧化磷酸化或糖酵解有關(guān)基因的均化表達量根據(jù)hOA進展分級制成熱圖,并對其均化表達量進行累計以直觀地比較。結(jié)果顯示,與軟骨細胞中氧化磷酸化相關(guān)基因的表達水平隨hOA進展升高(圖4C—D),而糖酵解相關(guān)基因表達水平隨著hOA進展而降低(圖4E—F)。

D:各級細胞氧化磷酸化相關(guān)基因均化表達累計值柱狀圖;E—F:各級細胞糖酵解相關(guān)基因均化表達水平熱圖及均化表達累計值柱狀圖。

A:氧化磷酸化相關(guān)基因沿擬時軌跡軸的表達模式;B:各級軟骨細胞處于S/G2/M期的百分比;C:各級細胞氧化磷酸化相關(guān)基因均化表達水平熱圖。

3 討論

關(guān)節(jié)軟骨是一種相對簡單的組織,僅由軟骨細胞組成。盡管如此,大量研究[10,19-20]揭示了軟骨細胞內(nèi)部異質(zhì)性的存在。scRNA-seq通過高效降維和擬時軌跡分析,可精確地對細胞內(nèi)部異質(zhì)性進行識別,在細胞譜系的建立和疾病發(fā)生發(fā)展機制的探索中發(fā)揮重要作用。有相關(guān)研究[10,21]利用scRNA-seq對hOA軟骨細胞內(nèi)部異質(zhì)性及其分子特征進行了探索,但對病情進展中基因表達水平及關(guān)鍵細胞生物學(xué)過程的變化進行分析。因此,本研究以骨關(guān)節(jié)炎連續(xù)性分級為基礎(chǔ),聯(lián)合單細胞測序擬時軌跡分析,挖掘并總結(jié)歸納出hOA病情進展中4類不同的生物學(xué)過程(圖2G),并對hOA病情進展中轉(zhuǎn)錄調(diào)控因子的變化趨勢進行了解析。此外,本研究還發(fā)現(xiàn),隨著hOA病情進展軟骨細胞代謝發(fā)生改變,表現(xiàn)為軟骨細胞氧化磷酸化水平加強,而細胞糖酵解程度降低。

骨關(guān)節(jié)炎是一種以進行性軟骨退化和關(guān)節(jié)炎癥為特征的疾病[22-23]。除炎癥反應(yīng)外,一氧化氮的生物合成、軟骨細胞分化、骨礦化、血管生成和入侵等生物學(xué)過程都與骨關(guān)節(jié)炎相關(guān)[22,24],本研究的分析結(jié)果也支持這一結(jié)論。與既往體內(nèi)動物實驗或體外細胞實驗不同的是,本研究基于骨關(guān)節(jié)炎分級進行了scRNA-seq數(shù)據(jù)再挖掘,不僅揭示了hOA發(fā)生發(fā)展過程中的相關(guān)生物學(xué)過程,還描繪出不同生物學(xué)過程發(fā)生的先后順序及交叉現(xiàn)象,為臨床開展針對性治療以阻止或減緩hOA病情進展提供了相對精準(zhǔn)的時間軸。

轉(zhuǎn)錄因子通過調(diào)控基因表達控制著細胞周期、細胞分化、細胞衰老和代謝平衡等重要生物學(xué)過程[25]。在hOA的發(fā)生發(fā)展過程中,轉(zhuǎn)錄因子調(diào)控因子的作用同樣至關(guān)重要。如RUNX2調(diào)控成骨細胞和軟骨細胞分化過程,在hOA發(fā)展過程中發(fā)揮重要作用[26],而SOX11過表達與軟骨損傷的惡化有關(guān)[27],本研究發(fā)現(xiàn)RUNX2或SOX11在關(guān)節(jié)軟骨損傷嚴(yán)重的后期表達顯著上調(diào)。IEZAKI等[16]發(fā)現(xiàn)轉(zhuǎn)錄因子ATF3在hOA組織中表達上調(diào),且ATF3敲除可減緩OA進展,但本研究結(jié)果提示ATF3的表達隨著OA進展而下調(diào)。產(chǎn)生這一矛盾的原因可能是組織與單細胞分析的異質(zhì)性。YANG等[28]對不同時期hOA的差異基因進行分析,發(fā)現(xiàn)ATF3隨疾病進展表達下調(diào),與本研究結(jié)論一致的。此外,除了已知的轉(zhuǎn)錄調(diào)控因子(如RUNX2、ATF13、SOX11和STAT1)[14-16,27,29],本研究發(fā)現(xiàn)了大量此前未知的在hOA進展中發(fā)生表達水平變化的轉(zhuǎn)錄調(diào)控因子,它們可能對疾病發(fā)展有調(diào)控作用。相較于傳統(tǒng)的基于對照組和hOA疾病組比較發(fā)現(xiàn)的差異轉(zhuǎn)錄因子,本研究利用患者軟骨細胞scRNA-seq數(shù)據(jù)發(fā)現(xiàn)的hOA不同發(fā)展階段的差異轉(zhuǎn)錄因子對hOA階段性靶向治療策略的研發(fā)有更加積極的意義。

除了生物學(xué)過程和轉(zhuǎn)錄因子,本研究還發(fā)現(xiàn)在hOA的發(fā)展過程中伴隨著細胞代謝的變化。細胞代謝不僅僅是細胞能量的來源,它還與細胞信號通路和表觀遺傳調(diào)控密切相關(guān)。此外,細胞代謝的變化也滿足了干細胞自我更新和分化等多方面的需求[30-31]。hOA的最終結(jié)果通常被認(rèn)為是軟骨祖細胞失去再生能力,逐漸朝向終端分化并失去自我修復(fù)能力[32]。這一現(xiàn)象是否與軟骨細胞的代謝改變相關(guān)目前尚不清楚。有研究[33]表明,hOA患者軟骨細胞中存在氮氧化物,這些氮氧化物會耗盡氧化磷酸化所需的物質(zhì)。此外,炎癥過程中產(chǎn)生的活性氧分子會影響呼吸鏈的運行,導(dǎo)致線粒體活性下降,從而影響軟骨細胞的氧化磷酸化水平[34]。還有研究[35]表明,hOA中軟骨細胞功能失調(diào)導(dǎo)致PFKFB3的表達下調(diào),而PFKFB3的表達水平則與糖酵解代謝強度呈正相關(guān)。以上證據(jù)表明,隨著hOA的發(fā)展,氧化磷酸化相關(guān)基因的表達水平上調(diào),而糖酵解相關(guān)基因的表達下調(diào)。另一項研究[36]發(fā)現(xiàn),線粒體DNA的單倍體G變異將導(dǎo)致細胞代謝從糖酵解轉(zhuǎn)向氧化磷酸化,從而增加了hOA的風(fēng)險。此外,hOA后期階段通常發(fā)生軟骨組織中血管生成,這會導(dǎo)致軟骨細胞所處環(huán)境中氧含量的升高,促使軟骨細胞從糖酵解代謝向氧化磷酸化代謝轉(zhuǎn)變。而相對于氧化磷酸化代謝,糖酵解代謝更有利于維持干細胞的自我更新能力[37]。本研究結(jié)果表明,hOA的發(fā)展伴隨著軟骨細胞中氧化磷酸化代謝的增強和軟骨祖細胞干性的逐漸喪失,這可能是hOA病程難以逆轉(zhuǎn)的關(guān)鍵原因之一,雖然這一假說尚需證據(jù)支持,但它表明代謝干預(yù)可能是阻遏或逆轉(zhuǎn)hOA進展的途徑,可作為后續(xù)臨床研究的方向。

綜上所述,基于hOA分級,本研究對軟骨細胞單細胞測序數(shù)據(jù)進行了再挖掘,探討了隨hOA病情進展的轉(zhuǎn)錄調(diào)控因子變化情況并歸納總結(jié)出hOA病情進展中的4種不同生物學(xué)過程演變,發(fā)現(xiàn)了疾病進展過程中伴隨的軟骨細胞代謝變化。本研究結(jié)果有望為臨床分階段的針對性治療hOA提供了基礎(chǔ),同時也為新療法的開發(fā)指明了方向。